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STUDY OF PHOTON CORRELATION TECHNIQUES 

FOR PROCESSING OF 
LASER VELOCIMETER SIGNALS 

By William T. Mayo, Jr. 

Summary 

The objective of this contract was to provide the theory 
and system design for a new type of photon-counting processor 
for low-level Dual Scatter laser velocimeter (LV) signals which 
would be capable of both the first-order measurements of mean- 
flow and turbulence-intensity and also the second order time 
statistics: cross-correlation, auto-correlation, and related 

spectra. 

This report provides a general Poisson process model for 
low-level LV signals and noise which is valid from the photon- 
resolved regime all the way to the limiting case of non-station 
ary Gaussian noise. Computer simulation algorithms and higher 
order statistical moment analysis of Poisson processes have 
been derived and applied to the analysis of photon correlation 
techniques. A Dual Correlate and Subtract frequency discrimina 
tor technique is postulated and analyzed. Expectation analysis 
indicates that the objective measurements are feasible. Error 
analysis for the mean-flow case indicates that practical 
transonic wind tunnel measurements are possible with 100-1000 
times less light than is required for burst-counter processors. 
A system design for a new high-speed photon processor for LV 
signals is provided. 
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INTRODUCTION 


The Problem 

Classical laser velocimeter (LV) electronic signal pro- 
cessing techniques are sometimes inadequate for detection of 
light scattered by small scattering particles which are required 
for following fluid motions. In other situations detection of 
larger scattering particles is difficult due to limited system 
sensitivity. Photon counting techniques offer improved system 
sensitivity by allowing velocity measurements to be made even 
when there are insufficient signal photons available to define 
the classical scattering signal. Such techniques are thus appli- 
cable when the presently used classical burst-counter and fre- 
quency tracker-techniques are not. 

The general objective of this contract was to provide the 
system design for a new type of photon-counting processor for 
low-level Dual Scatter LV signals which would be capable of both 
the first-order measurements of mean-flow and turbulence-inten- 
sity and also the second order time statistics: cross correla- 

tion, auto-correlation, and related spectra. This was to be 
accomplished by extending the preliminary feasibility analysis 
developed under a brief NASA Langley sponsored study* in early 
1974. In addition to theory, the system design would incorporate 
judgement based on experience in the experimental hardware devel- 
opment II] of a related, but simpler, photon-counting processing 
system designed and constructed for the U.S.A.F. Arnold Engi- 
neering Development Center to measure mean-flow velocities. 


The final report for that study (Contract NASl-13140) 
was informal and not disseminated. This report contains 
revised versions of all the necessary mathematics. 



Background 


Signal modeling .- Earlier modeling efforts have treated 
LV signals for which the noise could be considered as additive 
independent, stationary, and Gaussian [2, 3, 4, 5]. This is the 
limiting case of stationary Poisson shot noise which occurs for 
visible light photodetection when a steady light source such 
as a heterodyne reference beam [2] or high background light 
level [3] dominates the signal. In a recent simulation of low- 
level dual scatter signals, the accuracy of the noise model was 



Figure 1. Triply Stochastic Nature of Low-Level 
LV Signals; Turbulence, Bursts, and 
Photo electron Pulses. 
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extended by treating it as a nonstationary Gaussian process 
whose variance is proportional to the incident optical power. 

When we discuss "noise" in LV signal detection we are 
usually referring to the variation of the electronically 
detected signal with respect to a scaled version of the classi- 
cal optical (power) signal incident on the PMT. In a general 
analysis the classical LV signals are also random processes due 
to the random amplitudes and the arrival times of the signal 
bursts. Mayo [2] and Adrian [7] have treated these signals as 
a Poisson process for steady flow, and Durrani [4] and George 

[8] have treated them for the turbulent flow case as Gaussian in 
the limit of high particle number density. A new book by Snyder 

[9] treats generalized Poisson processes in great detail. 

Papoulis [10] provides a good introductory treatment. Snyder 
treats "doubly stochastic Poisson processes." These are inhomo- 
geneous (nonstationary) Poisson processes for which the rate 
function is a random process. Such a description is appropriate 
either for the classical LV signal bursts with the random turbu- 
lent flow affecting the rate of burst occurrence or for the 
single photoelectron pulses from the PMT with the random classi- 
cal bursts as the rate function. Clearly, when taken from the 
turbulence to the photo-electron pulses, a dual-scatter LV signal 
is a "triply stochastic" filtered Poisson process [11] . This 
three level nature of the signals is illustrated by Figure 1. 

Classical signals and burst counters . Presently accepted 
burst-counter and frequency-tracker LV processors were developed 
by analogy with wide-band frequency modulation (FM) and radar 
receivers. For FM and radar applications frequency detection 
(zero-crossing) circuits generally require about 10 db signal 
power to noise power ratio within the bandwidth of the system 
filters. The signals in such cases are continuous and the noise 

Modified version of noise model described in [6] . 
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is stationary and Gaussian, Comments by several speakers at 
the 1975 Minnesota IDA Symposium indicated that they had 
experimentally determined that their burst counter LV signal 
processors failed significantly when the signal power to nbise 
power ratio (during a burst) was less than 10 db. This 10 db 
condition occurs when there are approximately 10 or more photo- 
electron pulses per electronic response time Tj^. This response 
time is rise time or pulse width in the case of low-pass fil- 
ter. We have also defined this signal level as the lowest 
value of the "Gaussian" signal regime wherein the photomulti- 
plier current can be modeled as the classical signal plus non- 
stationary Gaussian noise; (although neither Snyder [9] or 
Papoulis [10] give any helpful rules as to when this asymptotic 
approximation is valid). For lower signal power the signals 
must be treated as Poisson. 

Photon resolved signals and photon counting .- A radically 
different approach to LV signal detection has been taken by 
Pike, Oliver, Jakeman, and others. Photon counting techniques 
were developed for use with low-level photon resolved signals. 
The summary results of several years of development of the 
single-clipping real-time photon correlator were described by 
Oliver and Jakeman in a recent book [12] . Dr. Pike described 
the application of photon-correlation to the processing of LV 
signals at the 1972 Purdue conference [13] . The presentation 
was apparently not received well by many attendees from the 
United States and little has been done in this country with the 
development of photon counting techniques until recently. 
Increased interest was shown by attendees of the 1974 Purdue 
Conference. 

One reason that the single-clipping correlator has been 
slow to acceptance in this country is that the original theory 
for its use was based on the assumption of many scatterers in 
the probe volume with the central limit theorem invoked to 
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render the statistics of the scattered electric fields Gaussian. 
Since this assumption is known not to predominate in many appli- 
cations of fringe-type LV systems in air, none of the first 
theory for a single-clipping correlator was directly applicable. 
Another problem with the existing commercial correlators for 
high-speed air flow besides speed (minimum time resolution of 
50 nsec) is the lack of any straightforward way to extend the 
concept to the determination of flow time statistics such as 
correlations and power spectra. The system proposed in this 
report addresses this latter deficiency as well as eliminating 
the problems of interpretation of single-clipping by using full 
multiplication. 

Several recent references provide additional valuable 
background information on photon correlation. Durrani and 
Created [14] provide a derivation of the expected value of a 
photon correlation from single particle (Poisson) signals. 

Birch et al. [15] have made experimental measurements in turbu- 
lent jet flows with skewed probability density functions. 

Abbiss et al. [16] also provide an analysis which shows that in 
some cases the Fourier transform of the correlogram may be 
interpreted as the probability density of the fluid velocity 
component. Durrani and Created [17] have investigated the use 
of some of the newer spectral estimation techniques which allow 
greater resolution from the limited number of data points in a 
typical correlogram. Finally, the reader should be aware that a 
new photon correlator instrument has been developed which was 
described by C. Fog [18] but the minimum time resolution inter- 
val is 160 nsec which is rather slow for high speed wind tunnel 
applications. (It was used for atmospheric studies.) 


We do not believe this will be a good approximation in 
many practical cases. This is discussed in a later section. 
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Burst rate/amplitude distributions .- The statistical dis- 
tribution of the classical burst amplitudes and the rate of 
occurrence versus amplitude are very significant in the charac- 
terization of any LV signal processor. It is generally accepted, 
for example, that the optimum rate of occurrence of bursts for 
a burst counter is less than the inverse burst duration (non- 
overlapping bursts.) It is also generally known that the error 
check circuits cause a burst counter to emphasize larger ampli- 
tude (good signal to noise ratio) signals. On the other hand, 
a photon counting processor must emphasize the lower amplitude 
signals in a distribution; the higher amplitude signals would 
produce only a single threshold crossing and otherwise be 
neglected by the system. It is therefore not possible to compare 
two different types of signal processors without knowing the sig- 
nal amplitude distributions and the processor behavior as a 
function of burst amplitude and other factors. Finally, in 
order to relate processor behavioi* to a specified particle size 
distribution, one must first relate the particle size distribu- 
tion to the burst amplitude distributions and then do all the 
other things already discussed. 

During the initial phases of our recent work for the USAF 
Arnold Center we addressed such questions as are suggested by 
the above statements both with theoretical models and with 
experimental measurements of burst rate/amplitude distributions 
for natural laboratory (unfiltered and unseeded) air. The reader 
is referred to the final report [1] for details. The following 
is one of the concluding paragraphs of that report: 

The statistical distribution of the amplitudes and 
rates of occurrence of classical bursts has been shown 
to be central in the problem of specifying or predicting 
the data rates and errors from any type of LV signal 
processor. Differential and cumulative rate/amplitude 
distributions have been formulated and analyzed theo- 
retically and have been measured experimentally for an 
argon backscatter LV system. The results indicate that, 
for the data obtained, the smaller aerosols contribute 
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more to the photon correlation accumulator than the 
larger ones. For the data measured, there would 
have been available less than 300 signals per second 
adequate in magnitude to produce burst counter data 
from scatterers larger than 0.7 ym in diameter 
while there would have been over 100,000 signals per 
second producing photon resolved signals from 0.2 - 
0.3 ym diameter particles. 

Scope 

In what follows we first develop general Poisson models 
for LV signals which include the non-stationary Poisson occur- 
rence of photo-electron pulses and the random amplitude effects 
of both the photomultiplier tube and the particle scattering 
cross sections. Formulas are provided for conditional and 
unconditional moments including mean, variance, auto-covariance, 
and higher moments (Appendix A). These formulas are for signal 
current, but they become valid for photon counting by a suit- 
able choice of the PMT output filter impulse response function. 

The next section evaluates the theoretical expressions 
for a specific Gaussian burst LV waveform model. These results 
are used to obtain the expected value of a photon correlation 
estimate. In addition a Dual Correlate and Subtract estimator 
which behaves as a statistical frequency discriminator is postu- 
lated and analyzed. The following section is devoted to statis- 
tical error analysis of the mean flow estimation technique using 
the Dual Correlate approach. The section after that shows that 
the statistical frequency discriminator may be applied to the 
estimation of turbulence correlations even though the time hist- 
ory of the velocity fluctuations is not available except as a 
noisy randomly sampled waveform. 

The results of the theoretical considerations and the exper- 
ience we have had previously with the AEDC [1] hardware study 
were utilized in a system design which is provided in Appendix 
D. Appendix C is a derivation needed in the section on variabil- 
ity errors. Appendix B provides the theory and an example pro- 
gram for correct Poisson simulation of low-level LV signals for 
evaluation of electronic processor models. 
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LIST OF SYMBOLS 


ENGLISH 


B 

Cii(t^,t2) 


f 

f 

m 

f( ) 


h 

h(t) 

H(o)) 

i 

i(t ) 
I 


J 

k 


m^(t) 

M 

.P 

M 

pq 

n(ti, tg) 

N 

N^ 


bandwidth 

autocovariance of i(t) 
electronic charge 
frequency 

mean Doppler frequency 

optical probe volume response function; pedastal 
function, peak amplitude normalized to unity 

random single-photoelectron charge gain 

Planck's constant 

impulse response of PMT and succeeding filters 
Fourier transform of h(t) 
counting index 

photocurrent at anode or succeeding filter output 
optical intensity 
counting index 
counting index 

time varying statistical mean current <i(t)> 

n, n, - n, n, or other function of n, 
k k-p k k-q k 

photon correlation sum at delay PAt 
accumulation of (pAT,qAx) for Dual Correlate Mode 
photoelectron count in interval (t^jtg) 
photon count in At interval about kAx 
number of time increments Ax (T = NAx) 

2 

number of fringes in the transmitter defined 1/e 
probe volume 


N. 


total number of AT intervals in advanced concept 
operation 


p integer delay number in photon correlation; 

largest delay number in dual correlate and subtract 

Pj^ constant optical background power 

F(t) optical power incident on photocathode 


q 


integer value of delay 
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R 

R(t) 

s») 

t 

"i 

T 

m 

u(t) 

U(t) 

U 

V(t) 

z(t) 


location of the nearest approach of scatterer to 
center of probe voliime 

constant rate of signal bursts 

instantaneous rate of signal bursts 

autocorrelation of X(t) =<X(t)X(t+x)> 

a generalized Poisson shot noise signal (Appendix A) 

time 

instants of photoelectron emission 

total data collection time 

inverse of mean Doppler frequency 1/fjjj 

time varying velocity component 

total velocity component 

mean velocity component 

vector velocity 
2 

l/e intensity radius at beam waist 
Poisson impulse process 


GREEK 

a 

3 

6(t) 

At 

At 

AT 


integration dummy variable; also l/e half-width 
of burst 

integration dummy variable 
Dirac delta or unit impulse function 
simulation time resolution interval 
photon processor counting interval 

accumulation interval for second level correlation 


e 

X(t) 



error 

time varying mean rate of photoelectron pulses 
X(t) = nP(t)/hv 

steady background photoelectron rate 

peak photoelectron rate (pedastal) from the jth 
scatterer 

optical wavelength 

mean signal photoelectron rate f <X.> 

product of photocathode quantum efficiency .and 
dynode collection efficiency 
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V 


p 


a^(t) 

T 

''h 
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(I) 

0 ), 

0) 


m 


optical frequency 

fractional turbulence intensity a /u * a /w 

u' u' m 

rms deviation of the radian Doppler frequencies 
from scatterers 

2 

time varying variance of i(t): <(i(t) - <i(t)>) > 
delay variable 

rise time or pulse width of low pass h(t) 

occurrence time for jth signal burst 

beam intersection angle; spherical declination 
angle from direction of incident light 

polar angle 

Poisson parameter 

Fourier transform variable or radian frequency 
random frequency of jth burst 
mean of random variable Wj (2‘Frf^) 


Special Notation 

oo 

<x> statistical expectation of x = / xp (x)dx 

-oo ^ 

OO 

f(t)*g(t) convolution: / f(a)g(t - a)da 

— OO 

^ denotes an estimate of a statistical average 


10 



STATISTICAL THEORY OF DUAL SCATTER SIGNALS 


Poisson Models 

T he Signal Current .- The signal current from a photomulti- 
plier tube (PMT) is modeled as inhomogeneous filtered Poisson 
random process (see Appendix A and also reference [19] ) given by 

oo 

i(t) = I eg.h(t - t.) (1) 

— OO 

where t^ = random time of the ith photoelectron 

e = electronic charge 

g. = random charge gain of PMT 

h(t) = impulse response of PMT/f liter system 

The system response h(t) is obtained as a convolution of the PMT 

impulse response h (t), the transmission line impulse response 

P 

h^(t), and the linear filter impulse response h^(t) 

h(t) = hp(t)*h^(t)*h^(t) (2) 

where the asterisk denotes the convolution integral: 

f(t)*g(t) = / f(a)g(t-a)da (3) 


The superposition assumes operation in the linear range of the 
PMT electron multiplier. The use of the function ^ (t) assumes 

ir 

that all single photo-electron pulses have the same shape except 
for amplitude. This neglects minor random shape variation. 

The quantity which relates i(t) to the classical optical 
power is the statistical mean rate X(t) of occurrence of the 
randomly occurring photoelectron pulses. Thus 

= nPm (4) 
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where 


n = product of cathode quantum efficiency and the 
dynode collection efficiency 

hv = Photon energy 

P(t) = Classical optical power, including background 

light and a constant Component for dark current. 

The effects of dark current are included by adding an 
equivalent power P^. Th^ model could be made more exact by add- 
ing a separate dark pulse summation with a separate distribution 
of amplitudes which are distributed somewhat differently than 
gjj^; but this distinction will not often be required in LV appli- 
cations. 

The previous material includes little which restricts it 
to LV signals. We now consider the form of X(t) which is also 
treated as a filtered Poisson process. 

Superposition of classical single burst signals .- Rigorous 
electromagnetic theory analysis of the scattered fields from 
more than one scatterer in the probe volume shows [20] mixing 
terms in P(t), the classical power incident on the PMT. However, 
in typical dual-scatter systems, the diffraction limited spot 
size of the collecting lens is much smaller than the probe vol- 
ume; conservation of energy arguments show that in such cases 
the number of scatterers in the probe volume may be much greater 
than unity with statistically negligible coherent mixing, regard- 
less of the quality of the collecting lens [2] . This is signif- 
icant even for LV systems which only trigger on isolated large 
signal bursts because we must also include in the model the 
effects of smaller scatterers which may exist at higher number 
density. We will take the position that at the PMT the classi- 
cal power P(t) is the superposition of the background light 
power and the power from individual scatterers without coherent 
mixing cross terms. This will be acceptable so long as the aver- 
age number of scatterers in one diffraction limited resolution 
cell of the receiver is less than unity. 
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A second consideration concerns the background light. 

Even when we neglect coherent mixing of signals, there are fluc- 
tuations in the classical background power. Bertolotti [21] 
provides a review of these effects. Broadband background sources 
can be largely suppressed by the use of narrowband spatial and 
wave-length filters, but not always adequately enough for meas- 
urements from small scatterers. If the background is modulated 
(for example fluorescent lights) the mean value signal is easily 
removed by electronic filters, but the non-statibnary noise is 
not. When the broadband background is "steady" there are actu- 
ally significant classical fluctuations at rates up to the opti- 
cal bandwidth. Bertolotti shows that when the optical filter 
bandwidth is much greater than the PMT electronic bandwidth, the 
photoelectron statistics behave as though the classical fluctua- 
tions did not exist (they are averaged out). Laser light scat- 
tered from windows is not broadband and may exhibit undesirable 
fluctuations. This background should be minimized, and its 
effects studied further. 


The random process X(t) .- With cognizance of the preceding 
discussion, we model P(t) as the summation of a constant back- 
ground P^ which includes broadband, laser, and dark current 
sources and an inhomogeneous filtered Poisson signal process. 


hvA ( t ) 

n 


P(t) 



oo 


+ I X f(t- 
d=-“ r 





hv 

n 


(5) 


where t. = occurrence time of jth scatterer reaching r., 

X. = random peak amplitude parameter, 

i3 

V. = vector velocity of the jth scatterer, 

r. = location of nearest approach of the scatterer 
^ trajectory to the center of the probe volume, and 

f(t,V,r) = normalized optical system response function. 

The notation in equation (5) explicitly shows that in general 

the shape of a burst (including signal period, signal envelope, 
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and pedestal) is a function of the scatterer vector velocity and 
trajectory location. The response also has random X. amplitude 
which depends on both trajectory and particle scattering cross 
section. In general the set of instants t. are independent 
Poisson random events whose instantaneous rate, R(t), is statis- 
tically correlated with the velocity vector. 

Discussion of the model .- Equation (5) is cast in a gen- 
eral form which obscures certain details with generality. First 
it assumes that the velocity of a scatterer remains constant 

while in the probe volume with a value V(x.,r.). The extended 

J 3 

theory of filtered Poisson processes is sufficiently general to 
encompass the fact that the functional form of the optical 
response function f(t,V,r) depends on two vector random varia- 
bles.* However, Snyder [9] assumes that the vector random 
parameters are independent. We are not certain at the present 
time what the statistical dependence of the rate function R(t) 
on the velocity V(t) may imply, but no serious consequences will 
result with low turbulence flow. 

Conditional Signal Statistics of the Photocurrent 

At times the models for the systems analysis problem may 
be simplified until analytical methods are applicable. In these 
cases the use of conditional statistics will usually simplify 
the analysis. Papoulis [10] discusses the use of conditional 
statistics at length. We utilize this technique at length in a 
later section. Basically for a multilevel random process the 
technique consists of assuming the higher level random processes 
are known and deterministic, evaluating conditional expectations 
assuming the higher level processes, then evaluating the expec- 
tation of the result with respect to the higher level processes. 
First we will consider statistics of i(t) assuming the classical 
optical signal X(t) is known. 

*Elementary shot noise theory is restricted to an impulse 
response function which is constant in shape. 
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Instan taneous mean, variance, autocovariance .- The result 
in Appendix A may be applied to determine the instantaneous mean, 
variance and auto-covariance of the signal. These are given in 
terms of the function A(t). The results are as follow: 

m^(t) = <i(t)> = e<gjj^>X(t )*h(t) 

= <(i(t) - <i(t)>)^> = e^<g^^ >X(t)*h^(t) 

^ii^^l’^2^ = <i(t^)i(t2>> - <i(t^)xi(t2)> 

oo 

= e^<g^^ > / X(a)h(tj^ - a)h(t2 - a)da 

CO 

where < > denotes statistical expectation and where the asterisk 

again denotes the convolution integral. These results include 

the specification that h (t), the impulse response of the PMT 

P 

anode, have unit weight, i.e., 

OO 

/ h (t)dt =1 (9) 

— 00 ^ 

in order to maintain conservation of charge. The functions 
h^(t) and h^(t) may include amplification or loss factors and 
need not have unity weight. 

Condit ional noise and SNR .- The concept of signal-to-noise 
ratio arose in communications theory when the "noise” was an 
additive stationary Gaussian random process totally character- 
ized by a mean, mean-square deviation (variance), and a power 
spectral density. The ratio of the peak or average signal to 
the rms noise was a useful measure. The preceding equations 
show the mean-square deviation (variance) to be an instantaneous 
time function which is related to the classical signal. Observa- 
tion of real LV signals on an oscilloscope display or computer 


( 6 ) 

(7) 

( 8 ) 


15 



simulations such as that shown in Figure 2 show that the concept 
of signal-to- noise-ratio is not an adequate figure of merit 
without careful specification. 

For an example SNR definition, we consider a low-pass PMT 
impulse response as a rectangular function: 

h(t) = ~ Rect (t/x. ) (10) 

^h “ 

where 

Rect (t) = l, lt|_<0.5 (11) 

= 0, It] >0.5 

If we now also assume that x, << T where T is the signal 

h m m 



Figure 2. Computer Simulation of LV Signal 

using Algorithms similar to Appendix 
B. (By J. F. Meyers NASA Langley.) 
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period of interest, then we could obtain an instantaneous SNR 
from equations (6) and (7) as 


SNR(t) = 


(t)/a^^ (t) 


<gi> 


<gi > 




<gi>V(t)/<g.^ X(t)/Tj^ 


( 12 ) 


2 2 * 
where the quantity <gj^> /"^^i ^ typically between 0.5 and 

1.0 with magnitude depending on the relative variance of the 

PMT single photoelectron pulse gain. For an ideal tube the 

quantity X(t)Tj^ would be the instantaneous SNR. This is not 

useful since it is a time function instead of a number. 


As an alternative, we may take the local time average of 
the SNR given by (12) over a single cycle near the peak of the 
pedastal. This would give, for an ideal PMT, 


S>®AVpeak ” 

where X. is the peak value of the pedastal of the jth signal 
burst, if we assume sparse non-overlapping bursts. We observe 
that equations (6) through (8) are valid when h(t) is a bandpass 
function, but (13) is then meaningless unless we redefine 
for a bandpass h(t). Also we note that this definition would 
be necessary for meaningful use with a burst-counter processor, 
since it is the bandpass filtered AC signal to wide-band noise 
power that is significant in that case. 

Signal Regimes .- The idealized quantity ®NR^^ given 

in equation (13) as is at least a useful quantity in defin- 

ing a classification system of signal regimes for a low-pass 


♦ 


Typical rms values are 0.707 or greater ([22], page 66). 
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filtered signal. The following definitions of a photon-resolved, 
a photon-limited, and a Gaussian signal regime have been some- 
what arbitrarily* identified: 

The signal is photon resolved if << 1. In this case 

the probability of two or more photoelectron events occurring 
within the response time is small. Its appearance is that 
of individual pulses which vary in height due to the randomness 
of g^. Photon counting methods are appropriate. The condi- 
tional mean value of i(t) is still proportional to X(t), but 
there is no visible resemblance to X(t). This condition is 
illustrated by the extreme right hand portion of Figure 2. 

For XXj^ >> 1, the signal i(t) is asymptotically a nonsta- 
tionary Gaussian Process. In this case the first and higher 
order probability density functions for i(t) at any set of 
instants (tj^,t 2 ,...) may be determined immediately by plugging 
the mean, variance and auto-covariance from the preceding equa- 
tions into well known Gaussian formulas. Under these same con- 
ditions the signal display appears to the eye as a classical 
signal m^(t) plus Gaussian noise. This condition is approached 
by peak of the trace in Figure 2. The major difference between 
this case and that of classical communications theory problems 
is that the a value for the noise is signal (time) and system 
dependent. Usually, signals in the Gaussian regime are suitable 
for processing by classical methods (burst counter and/or 
tracker ) . 

The photon limited regime is that for which Xx^^ is within 
an order of magnitude of unity. No mathematical simplifications 
are possible. Visually the signal appears as shown in all but 
the lowest portions of Figure 2. The upper limits of photon 


See Papoulis [10] page 571. No specific limit on the 
magnitude of Xx^^ is given. We define the photon limited regime 
as 0.1 < XXjj < 10. 
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counting techniques and the lower limits of conventional tech- 
niques both fall in this range. 

As we have illustrated, the signal classification may 
apply to different portions of the same waveform. We may also 
use the classification with respect to the peak pedastal value 
to classify signal bursts. Under this type of classifi- 
cation, Figure 2 illustrates a photon-limited burst whose peak 
average SNR is less than 10. Additional bandpass filtering 
would increase and place the central portion of the burst in 
Figure 2 in the Gaussian regime. This would not be possible 
with a significantly weaker signal. 

Unconditional Statistics 

L ong time mean, variance, autocovariance .- Equations (6) - 
(8) include the assumed deterministic classical signal X(t) which 
is proportional to instantaneous optical power. When we wish 
later to evaluate the long-time average result which accumulates 
during a photon counting experiment, it will be necessary to 
treat X(t) as an ergodic random process with long-time average 
equal to the unconditioned statistical mean; 

<X(t)> = X (14) 

We also make use of the autocorrelation of X(t): 

<X(t) X(t + T)> = Rj^(T) (15) 

Now from equation (6) taking the expectation with respect to 
X(t) gives the average current as 

oo 

<i> = e<g.xX> / h(t)dt (16) 

^ —00 

Id 



where the integral is unity unless h(t) includes preamplifi- 
cation or attenuation external to the PUT. In order to deter- 

2 

mine the long time variance cr. we do not get the correct 

^ 2 

answer by taking the expectation of a. (t) given by equation ' 

^ 2 - 

(7). Rather, one determines the conditional value of <i (t)> 

by adding the square of (6) to (7). The expectation with 

respect to X(t) follows; finally, the square of equation (16) 

2 

is subtracted from the unconditional expectation of i (t). When 
all these steps are completed, and similar ones for the uncondi- 
tional autocovariance, we obtain: 

oo oo 

= e^[<g^>^ / R^(a)fj^(a)da + <s^><\> / h^(a)da (17) 


- <gj^>^<X>^ (/ h(a)da)^] 

— OO 


Cii(x) = e^[<g^>^R^(T)*fj^(x) + <g^^xX>fj^(x) ( 18 ) 

OO 

- <g.>^<x>^(/ h(a)da)^] 

— OO 


where 


fjj(x) = / h(a)h(a + x)da (19) 

— OO 

The second term in the expression for C^^(x) vanishes for x 
greater than the impulse response time for the PMT and filter 
combination; the last term is the square of the mean; the first 
term is the correlation of X(t) smoothed by the correlation of 
h(t) with itself. 

Ideal photon correlation .- An idealized photon correlator 
counts all photoelectron emission events during successive uni- 
formly spaced clock periods of duration Ax. The number sequence 
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which results is algebraically manipulated to yield the 


summation of terms 


In evaluating the expected value 
of the result of acciimulating such a sum we encounter the need 
to evaluate the quantities <nj^> and + ^k+p^‘ expec- 

tations may be evaluated using equations (16) - (19) by assuming 


h(t) = Rect(t/AT) (20) 

Si = 1/e 

where Rect(t) was defined in equation (11) and 

/h(t)dt = At (21) 


With these assumptions, i(t) is equal to the number of photo- 
electron events in the interval (t- At/2, t + AT/2) and the 
formulas reduce to 


<i> = = <X>At 

oo 

var = <n^ = At / Rj^(a) A(^)da 

— -OO 

+ <X>At - <X>^At^ 


( 22 ) 

(23) 


..2 


<Vk+p> = C. .(pat) h- <i> 


= ^'^Rx('^>*^('^/^'^>It=pAt,p?40 


(24) 


where the correlation integral of equation (19) produces a tri- 
angular function, i.e. 


/ Rect (^) Rect(^^)da = AtA(t/At) (25) 
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where 


A(t/At) = 1 - |t1/At , It I < At (26) 

= 0 , I T 1 > At 

We observe that the generally accepted result that the photon 
correlation is shaped like the correlation of the classical 
signal is true subject to the jump discontinuity at zero delay 
and subject to the triangular weighting function which behaves 
as a low-pass filter with respect to the details of 

When At is much smaller than a characteristic signal per- 
iod, then equations (23) and (24) simplify: 

^Vk+p^ = R^(PAt), ^ f 0 (27) 

= At^ R^(0) + <X>At, p=0 

var n^ = = <X>At + Ax^ (<X^> - <X>^) (28) 


i 


e. , 


2 

var nj^ = + Ax var X (29) 

This last result, which we obtain as a special case, has also 
been given by Bertolotti [21] . It provides a way to measure the 
variance of the classical signal even with photon resolved sig- 
nals if a long sequence of nj^ values are available. 

PHOTON COUNTING PROCESSORS FOR MEAN 
FLOW AND TURBULENCE INTENSITY 

In this section we provide an idealized theoretical basis 
for the use of photon correlation and a new type of photon 
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counting technique, "dual-correlate-and-subtract , " and show how 
they may be used for mean-flow and turbulence intensity measure- 
ments. Extensions to higher order statistics are discussed in a 
later section. Error analysis of the processors described in 
this section is presented in the next section. 

Specific Signal Model : 

The form of equation (5) is quite general. In most prac- 
tical systems, the function f(t,V. ,r.) which describes the opti- 
cal response with respect to particle position and velocity is 
complicated when the effects of limiting pinhole apertures and 
variable duration due to high turbulence are included. For the 
present we assume a simplified low turbulence model which assumes 
a burst with perfect contrast and constant shape: 


X(t) = + I Xjf(t - Tj)tl + COS 0)j(t-Tj)l (30) 


where X . = random burst pedastal amplitude 

Xj^ = constant background rate, 

f(t) = low-pass pulse waveform with peak equal unity, 

X. = occurrence time for jth burst, 

(D . = radian frequency proportional to one velocity 
^ component of the jth particle 


We may write equation (30) as the sum of a constant, X^^, a low- 
pass process, X^p(t), (the pedastals), and a bandpass process 
Xbp(t)- Then 




(31) 


We assume that there are several fringes in the probe volume so 
that the spectra of X^^(t) and Xj^p(t) are non-overlapping. Thus 
Xbp(t) is a zero-mean process, and X^^p(t) and X^p(t) are uncor- 
related. We obtain 
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( 32 ) 


B,(t) = <X(t)>2 H. 0„p(T) * C^j,p(T) 

We now use the low turbulence assumption and further assume 
that the scatterers are uniformly dispersed in space so that, 
R(t), the rate of burst arrivals is a constant R. The results 
of Appendix A can be applied to derive expressions for the 
three terms in equation (32), the result is 

00 

<X(t)> = + R<x.> / f(t)dt (33) 

— no 


Cj^^p(T) = R<X^> / f(t) f(t+T)dt 

oo 

Cj^bp(’’^) ~ ^ <COS WT> / f (t )f ( t +.T )dt 


(34) 

(35) 


where the expectation of the <cos wx> term is with respect to 
the random variable w . . The derivation requires that we expand 

tJ 

the product of cosines with the sum and difference formula and 
approximate the integral of the product of a low-pass term and 
a bandpass term as zero. 


We now assume that the turbulence is Gaussian with mean 
radian frequency and rms deviation o^. Then by direct appli- 
cation of the definition of expectation we obtain* 


<COS 0)T> = / 


-(to-w )^/2a^ 
m ^ (A) 


cos 0)T do) 


«> / 2tt a 


0) 


(36) 


= e 


2 2,_ 
-a T /2 
w ' 


cos (I) T 

m 


*In general <cos oot> is simply related to the statistical 
characteristic function for the random variable oi . since this 
function is defined as a Fourier transform of the'^probability 
density function [10] . 
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A simplified expression for the autocorrelation of the classical 
signal is thus 


R 


xC-r) = + R<^j> / f(t)dt)' 


(37) 


1 -^ 2/2 


+ R<x 2>(1 + i e “ 

J ^ 


cos (*) t)J f(t)f(t+T)dt 

— 00 


This result shows that regardless of the shape of the envelope 
function f(t), the autocorrelation function has a consinusoidal 
variation at the mean signal frequency and a Gaussian envelope 
decay factor due to turbulence intensity. 

If we now assume that the classical bursts are Gaussian 
shaped (TEM^^ beams without aperture effects), then we have 


f(t) = e 


. 2 , 2 

-t /a 


(38) 


where a is the 1/e half width of the envelope and obtain 


oo 

/ f(t)dt = /Tr a 

— 00 


OO 

/ f(t)f(t + T)dt 

— OO 



/2a2 


(39) 


The final simplified expressions for the first two moments of 
X(t) are 

<X> - <X(t)> = A, + R<A.> /Ir a (40) 

^ 0 

^ 2 ^ 2^2 2 2 

R^(t) = <A>2 + R<x2>y|'ot(l + I e “ cos o)jjjT)e“'*^ /2ct ^ 43 ^^ 

where low percent turbulence has been assumed and 
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= background count rate, 

R = scatterer arrival rate, 

Xj = random pedastal height from jth scatterer, 

a = 1/e half burst duration, 

= rms deviation of radian frequency due to turbulence, 

m « mean radian frequency (a /w_ « 1), and 
m 0) m 

T = delay variable of autocorrelation. 

The shapes of typical correlation functions for zero tur- 
bulence and 10% turbulence are illustrated in Figure 3. 


Idealized Photon Correlation of LV Signals 

The number is the number of photo-electron emission 
events in the interval which extends ±At/2 from the 
instant kAx. An idealized photon correlator produces and 
sums delayed products from the uniformly spaced sequence 
We assume the total number of products accumulated is N. The 

y\ 

accumulator produces a sum at the delay value pAx given by 



Ji 


(42) 


The ideal photon correlator would simultaneously accumulate Nnj^ 
defined by 





(43) 


The unconditional expected value of these sums is obtained from 
equations (22) and (27) after interchanging expectation and sum- 
mation as 


<Mp> 


N<n, n, . 
k k+p 


> 


= NAx^Rj^(pAx), p 0 
= Nl<X>Ax + Ax^Rj^(O)], p=0 


(44) 
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<Niij^> = N<X>At 


(45) 


By combining equations (40) - (45), we may obtain an esti- 
mator for the autocorrelation of X(t), i.e.. 


N 


ft (pAx) = — I n n . , p f 0 
^ Ax^N Ax'^N k=l ^ ^ P 


(46) 


^x(PAx) “2^ "“p 

Ax N ^ 


(M^ - Nnj^) , P = 0 


N 


2 i “k'"“k 
Ax N k=l ^ ^ 


I P = 0 


This estimator includes the zero delay value, which is usually 
omitted, by making use of the separate mean count computation. 
The same mean count estimate nj^ may be used to estimate the long 
delay level 

<nj^/Ax>2 = <X(t)>^ = Rj^(~) (47) 

Interpretation of an autocorrelation estimate computed 
according to equation (46) involves, first, the use of an ana- 
lytical model such as equation (41); second, a parameter extrac- 
tion procedure, such as a mean square error minimizing curve fit 
algorithm; third, a correction for any statistical bias errors; 
and fourth, a variability error criterion which assures that 
sufficient data is accumulated. We have provided a procedure 
and an example model for the first step. The bias and varia- 
bility errors are discussed in a later section. We have not 
considered the optimization of the second step although one 
method is discussed below. Some literature [15,16,171 is begin- 
ning to appear, but further effort is needed. 
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At the Denmark LDA conference in the summer of 1975, 

Abiss et al. [16] described a "new" interpretation of a photon 
correlellogram as the Fourier transform of the velocity proba- 
bility density function. This "breakthrough" allows the meas- 
urement of the probability density of the velocity field by 
digital Fourier transform of the correlation results under the 
condition that the pedastal and fringe envelope correlations 
are nearly constant over the range of delays for which the sin- 
usoidal correlation is appreciable. This result is equivalent* 
to the interpretation of spectrum analyser displays (connected 
directly to the photo-detector) years ago as the probability 
density of the velocity field. The restriction is equivalent 
to requiring many fringes in the probe volume so that the transit- 
time spectral broadening is small. Our equation (35) shows the 
relationship between the probability density function of the 
velocity (frequency) samples and the correlation. The expres- 
sion <cos o)T> is the cosine Fourier transform of the probability 
density for coj . When there are many fringes in the probe volume 
the envelope correlation term is a broad pulse which is essen- 
tially constant over the extent of <cos (jjt>. In that case, the 
transform of the correlation is actually inverting the statis- 
tical characteristic function to produce a scaled probability 
density. 

The unfortunate truth remains that in most practical LV 
problems, the probe volume and spatial frequency of the fringes 
must be small with the result that we rarely have the luxury 
of having many fringes in the probe volume. For the realistic 
case, then, we find the spectrum to be the convolution of the 
velocity probability density function and the transform of the 
envelope correlation, and care is required in interpretation. 


At least in theory: the photon correlator is much more 

efficient at low signal levels than any swept frequency spectrum 
analyser could be. 
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A Photon Counting Frequency Discriminator 

This section describes the basis for a statistical mean 
frequency discriminator for photon resolved LV signals. We 
begin with a historical description of the motivating logic 
which may offer insight to others for more advanced development. 

Autocorrelation Frequency Discrimination .- It is well 
known that the autocorrelation function of a zero-mean narrow- 
band random process is cosinusoidal with frequency equal to 
the mean frequency of the narrow-band process and envelope 
which is the autocorrelation of the amplitude and phase envelope 
of the process. For example, if the two-sided power S (f) of 

A. 

a random process x(t) is 

+ (48) 

Where A(f) is an even pulse-shaped function and f^^^ is a center 
frequency much greater than the frequency width of A(f), then 

Where A(f) is the Fourier transform of R (t). Now we note that 
for small perturbations of the argument SuTf^^j about the point 

T = 3/4 f^ (50) 

the value of R (t) varies in proportion to the perturbation of 
either t or f^^^. Thus, if x is selected to satisfy equation 
(50), then the statistical autocorrelation function would serve 
instantaneously as a frequency discriminator lor small deviations 
- 

This is significant in later sections even though the 
statistical autocorrelation function cannot be experimentally 
observed. 
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of fjjj. If the random process were ergodic, then the same 
results would apply to time-average autocorrelations when the 
averaging time was long compared with the coherence time of the 
process. These principles may be extended to a random process 
which is the sum of a narrow-band process and a low-pass pro- 
cess by first filtering the process with a high-pass filter to 
remove the non-zero mean low-pass process. 

Appli cation to photon resolved signals .- In the case of 
photon-resolved LV signals, the classical signal random process 
is not necessarily recoverable from the sparse single photon 
events, but we have already shown that the autocorrelation 
function may be approximated by the expected value of a photon 
correlation operation on the photon events. The original con- 
cept was thus to devise a high-pass digital filter which would 
be applied to the count sequence tnj^} to remove the effects of 
background light and low-pass pedastal from the statistical 
average while leaving the bandpass information signal 

information. This filter would be followed by a digital corre- 
lation (delayed product summation) at the delay value given by 
equation (50). Since the digital electronics had to be very 
fast, only simple digital filters could be considered. The 
simplest one that we thought of was to delay the sequence 
by one-half period of the signal and subtract. This operation 
cancels the low-frequency portions of the expected value of 

and adds the sinusoidal portions in phase. Such a delay and 
subtract filter would produce 


m. 


= n, - n 


k-q 


( 51 ) 


Where the counting interval At must be adjusted to satisfy 


with q an even integer. After this was done one would accumu- 
late lag products of mj^ with delay pAx equal to three quarters 
of a signal period p = 3q/2. This approach leads one to form 
the summation given by 


/V 


M 



N 

y - n, n, + 2n, n, - n, n, 

k-p k-q k k-p k k-p-q 


(53) 


The previous results for photon correlation are applicable for 
the determination of the expected value of the quantity 
The details are omitted here and included for the simpler dual- 
correlate approach described below, but Figure 4 is a simplified 
discriminator characteristic for the expected value of the 
result. It is provided for comparison with Figure 5 which is 
described in detail below. 

After all of the above reasoning, it occurred to us that 
another approach (we thought) consisted of performing the delay 
and subtract filtering operation on the correlation estimate 
after it was made rather than on the high-speed signal sequence 
{n^}. Reference to Figure 3 shows that a one-half-period shift 
and subtract of the typical LV correlation will approximately 
cancel the low-pass portion of the autocorrelation. The results 
of this approach provide an approximate discriminator response 
as shown in Figure 5 with less high-speed data arithmetic 
required. This "dual correlate and subtract" approach also 
leads to larger frequency range for the discriminator function, 
and thus has been chosen for development. Hindsight has now 
shown us that by algebraic rearrangement this approach is most 
economically implemented with high-speed delay and subtract 
filter prior to multiplication. This technique may be shown 
to be just another implementation of the original filter and 
correlate concept. 
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Dual correlate and subtract .- We will let the quantity 
M be defined as follows 

pq 



(54) 


where 


”k * "k^“k-p - ”k-q> 


(55) 


It is straight forward to show that, except for a few end 

terms, M is mathematically identical to the quantity M -M 

where M and M are defined by equation (42). It is for this 
P Q. 

reason that we will lable the approach as the "Dual Correlate 
and Subtract" technique even though the delay values are nega- 
tive instead of positive.* We now demonstrate that the expected 

ys 

value of behaves as a frequency discriminator as illustrated 

pq 

in Figure 5 under conditions which we will identify. The 
adjustment of the system clock period Ax leads to a null in the 
expected accumulator value. Measurement of Ax provides a direct 
measure of the mean signal frequency as we shall now show with 
our simplified signal models. 

From equation (44) we obtain the expected value of the 
estimator as 


= NAx^[R, (pAx) - R, (qAx)] (56) 

pq A A 

The complete expression is obtained for our simplified signal 


♦ 

The negative delay implementation was more suitable for 
the hardware design. 
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model by using equation (41) in its entirety. Here we assume 
that pAt and qAx are both small compared with a so that the 
pedastal terms cancel as well as the steady term. This leaves 

<Mpq> = NAt^R<x|> I a • (57) 

• e cos 0 ) dAt - e cos w„qAT . 

m m 


Now we require that 


p = 3q (58) 

and observe the behavior of equation (57) near the values of 
At specified by letting qAx be one quarter of the signal per- 
iod where both cosine terms will vanish. 

qAx = 2 tt/4Ujjj = T^/4 (59) 

The shape of the term in braces is plotted in Figure 5 
under the assumption of many fringes in the probe volume 
(a large) and low turbulence (a^ small). Thus the quantity 
plotted is simply 


[cos 3qAx 0 )^ - cos qAx to^] (60) 

Figure 5 illustrates the expected behavior of the accumulator 
sum for changes in the mean signal frequency If the sys- 

tem clock frequency is changed to change Ax, then the response 

2 

is the product of Ax and the curve shown in the figure. The 
shape of the curve is affected but the zero crossing locations 
are not . 
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Selection of delay constants .- The theory does not 
uniquely specify the relationship of At (the above system clock 
period) to the signal period because q is not specified. 

For a given value of signal period, the largest possible 

value of At for an acceptable null is when q = 1 and 
At = Tjjj/4. This produces the least variability error (as we 
will show) and the most bias error due to time smear (the tri- 
angular weighting function in equation (24). The value q = 1 
also allows the highest signal frequencies to be measured for 
a given maximum system clock frequency. The bias errors can 
be reduced at the expense of increased variability error and 
reduced maximum signal frequency by using q = 2, At = T^/8. 

In order to facilitate experimental research, our design for a 
research instrument allows selection of p and q over a wide 
range. 


STATISTICAL ERROR ANALYSIS 
Introduction 

The two principle types of error which arise in statis- 
tical measurements are bias error and variability error. Bias 
error is a term which refers to the difference between the 
statistical expectation of the measurement system output and 
the desired average value being measured. The variability 
error is the rms value of the random deviation of a specific 
experimental result from the statistically expected value. For 
ergodic random processes, the variability error converges to 
zero in the limit of infinitely long data collection time; but 
it converges to an acceptable level (which must be defined in 
the measurement objectives) within a finite measurement time. 
The bias error cannot be removed by further averaging but it 
can often be removed by analytical compensation or by experi- 
mental calibration when it is small compared with the desired 
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quantity. In general, analysis of both types of errors is 
required in any statistical measurement. 

There are three different methods for evaluation of the 
statistical errors of a measurement system; these are: analysis, 

simulation, and experiment. Appendix A provides theory for the 
higher order moments of inhomogeneous Poisson processes which 
we have expanded for application to the demonstration of con- 
cepts in the previous sections and analysis of errors in this 
section. The study of the theory of Poisson processes has also 
provided the concepts necessary for digital simulation of the 
LV signals and their detection by photon counting systems. 

The concepts and a FORTRAN computer program which we have 
developed for this simulation are provided in Appendix B. The 
simulation provides a method of investigating such nonlinear 
effects as processor dead time and counter saturation which we 
have not yet been able to do with analysis. 

Bias Errors 

The sources of bias errors which have been studied in 
the measurement of mean flow are as follows: fringe number, 

turbulence intensity, time smear, and dead time. 

Fringe number and turbulence intensity .- When the number 
of fringes N in the probe volume is small, subtraction* of 
Rj^(t = 3T^/4) - is slightly negative instead of 

zero. The error becomes significant only for small values of 
Nj. The signal frequency estimate using the zero criterion 
will be too small. We have computed the error in percent in 
the following manner. 

Let Rj^(3qAx) - R^(qAx) 4 tJ^(qAx). For large nvimber of 

- 

This error pertains to the dual correlate and subtract 
technique. 
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T 

fringes (N^~) then i/^CqAx) = 0 when At = where is the 
mean signal period. When N«», then 


= Ae = i|;(At = Tj^/4) j< 0 (61) 

The correct value of qAr for 4/(qAx) = 0 is qAx = » (T /4) 

and the fractional error e in accepting ij^(qAx) = 0 as indica- 
tion that qAx = T^/4 is approximately 


4tj^(V4) 

T li’ (T /4) 
m^ ^ m' ' 


(62) 


We have evaluated e for various levels of p = rms turbulence/ 
mean velocity and for = number of static optical fringes in 
1/e^ signal width. 

The result of the parametric computations are presented 
in two forms in Figures 6 and 7. From the results we conclude 
that the bias error is not very sensitive to turbulence 
Intensity at the levels shown. A first order correction inde- 
pendent. of turbulence is thus possible. A second order correc- 
tion is possible with only very imprecise estimates of turbu- 
lence intensity. 

Time smear error and correction techniques .- The photons 
represented by the count nj^ at time kAx were in fact smeared 
in occurrence time over the range [(k-l/2)Ax, (k + l/2)Ax]. 

The effect of this is usually neglected by assuming Ax is small 
compared with any significant variations in the classical sig- 
nal. This is never true in the most difficult experimental 
cases where the electronic speed limitations force Ax to be an 
appreciable function of the signal period T^^^. The result pro- 
vided in equation (24) for photon correlation included the time 
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FRACTIONAL TURBULENCE (p ) 


Figure 6. Turbulence Intensity Mean Bias 
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smear effect in the triangular function A which convolves with 
the autocorrelation of ' the classical signal. After- equation 
(26) the triangular function was treated as a delta function, 
i.e., in the limit of small Ax, : 

^ A(t/At) ^ 6(x) (63) 

and equation (27) results from the convolution. 

Now, if we do not make the limit assumption,* we note 
that the effects of the triangle function are easily displayed 
in the frequency domain: From Bracewell [24] we have the rela- 

tionship . 

oo 

/ A(T/Ax)e“^^'^^ = sinc^ (fAx) (64) 

OO 

where 

/j,* ^ _ sin irfAx /nc\ 

sine (fAx) = (65) 

This Fourier transform relationship is illustrated in Figure 8. 

The convolution theorem of Fourier transform theory assures us 

that in the frequency domain, the effect of the convolution in 

equation (24) is a product. In other words, the frequency 

spectrum associated with the signal correlation by Fourier 

transform is attenuated by a low-pass filter whose form is 

sine (fAx). This function is plotted in Figure 8c. As the 

figure shows, there is little attenuation when f^Ax = Ax/T =1/16, 

m ' m 

i.e., when the clock period is 1/16 of the signal period. For 

ijc ' 2 

The effect of time discretization is a sine (f)type of 
low-pass filter. This is similar to a result in our previous 
work of determing turbulence power spectra from randomly tried 
samples [23] . 
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Signal Part of 


Weighting 

Functions 



(a) Triangular Weighting Functions corresponding to choices 
for At relative to mean signal period. (illustrated at 
delay locations for Dual Correlate Approach) 
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sine X 
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(b) Fourier transform relationship of A and sine function 


sinc^ (f At) 
m 



At/T = fAT 
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(c) Low Pass Filter effect of At selection. 


Figure 8. Time Smear Effects of Finite At. 



very low turbulence and a large number of fringes N^, the 

signal spectrum would be small in width and At could be made 

as large as T /2.* However, for values of At > T /8 we see 

that a correction of the velocity probability density function 

could be useful for all but very small turbulence levels. 

Such a correction would be effected by multiplying the trans- 

2 

form of the correlation estimate by 1/sinc (ATf) prior to 
final data interpretation. 

Dead time effects .- No physically realizable photo- 
detector and electronic counter combination can be constructed 
without some dead time; i.e., a period of time following a 
threshold crossing by the analog photodetector voltage wave- 
form during which no additional crossing events will be 
counted. This dead time is typically 10 nsec for commercial 
photon discriminator circuits at the time of this writing. 

There is no fundamental reason why this cannot be reduced to 
less than Snsec with the fastest photomultiplier tubes and 
counting circuits now available. We will also distinguish two 
other types of counting dead time. The first of these is 
"pulse pile up" in which the photodetector analog waveform 
remains above the threshold level due to there being more than 
one photoelectron event within the pulse response time of the 
photo detector. There is also a brief interval during the 
periodic counting interval during which the counter is being 
reset and is not available . This is true even when two counters 
are used alternately, since a small guard interval is required 
to keep both counters from counting the same event. Practical 
dead intervals t^ may be on the order of 2nsec or less. 


3|C 

This limit could be actually exceeded because the 
sampling theorem really refers to bandwidth, not maximum fre- 
quency; however, the attenuation would then be severe. 
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The counter dead interval will be considered first 
since it is the simplest. To first order it can be neglected, 
but correction in the formulas is simple. It is only neces- 
sary to replace At by At-t^ in the theoretical formulas where 
the At refers to the counting width. The amount of delay 
remains equal to some integer multiple of At. As a simple 
example the mean value of becomes <X(t)>(AT-T^) instead of 
<X(t)>AT. 

The effects of pulse pile up are least known for the 
photon-limited signal cases. For higher photon rates the PMT 
analog waveform will generally remain above the photon counting 
threshold and so produce no counting effects. The pulse pile 
up effects can be studied with the present signal simulation 
program (see Appendix B) when a more complete processor simu- 
lation subroutine is completed in the future. At present the 
photon processor simulation is idealized so that discretized 
photoelectron event times are used directly without synthesis 
of the PMT anode waveform and a threshold crossing circuit 
model. 

The effects of discriminator dead-time effects have been 
presented analytically by Jakeman* for the case of Gaussian 
optical electric field statistics. We have not yet extended 
the theory to the single-particle LV signal situation (non- 
Gaussian field statistics) and can offer no improvement analyti- 
cally here over Jakeman 's results. However, we dispute Jakeman's 
conclusion* that in practice dead time errors are not a serious 
limitation of photon correlation. The seriousness of the 
effect is demonstrated below using the computer simulation pro- 
gram provided in Appendix B. 

We have simulated a wind tunnel instrumentation 
example with the constants given in Table 1. The results are 

*Page 116, [12]. 
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Theory 

Symbol 

Program 

Name 

Value 

Description 

T 

TTOT 

10 ^sec 

Total time limit 


MNTOT 

3x10^ 

Photoelectron limit 

Enjj. 

KMAX 

480 

Number of At ' s per burst 

<Xj> 

AEAC 

10^/sec 

Mean peak pedastal rate 

HIGH 

10 

Max Xj = High*<Xj> 


LOW 

0.1 

Min X . = Low *<X . > 

J J 

Xj = AFAC (not random) 


IR2 

4 

1/R 

TB 

5 X 10"® 

Mean time between bursts 


THEORY 

TRUE 

Logical: Selects ideal 

processor 


MAXC 

10® 

Maximum count of processor 
(high to avoid limitation 
here) 

ct/2 

D 

8 X 10”^sec 

1/e® half width of bursts 


CON 

3 

D CON = Total burst width 


AME 

1.0 

Burst modulation index 

fn, 

m 

FO 

6.25 X 10®Hz 

Signal frequency 

Ar/At 

ITAU 

4 

Ratio of processor resolu- 
tion to simulation resolu- 
tion 


IP 

20 

Maximum correlation delay 


WAVE 

TRUE 

Logical: True=Bursts present 

Xb = o 

CONST 

FALSE 

Background : True=X^ present 

Xd/At 

DEAD 

0,1, 2, 3, 4 

Dead time in At units 

At 

DT 

5 nsec 

Simulation resolution 
interval 

At 

DT+ITAU 

20 nsec 

Processor resolution 
interval 


ONE 

FALSE 

Logical : Time=Dual Correlate 

TABLE 

1 . PARAMETER 

SELECTIONS FOR 

DEAD TIME SIMULATION. 
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plotted in Figure 9. For comparison, the figure also shows 
the theoretical expected value of the accumulator sums as 
predicted by the following equations. We use the values in 
Table 1 with equations (40), (41) and (44) to obtain the 
mean level <X(t)> as 

g 

<X> = 2.005 X 10 photoelectrons/sec (66) 


— 7 2 

<M >= 80.4 + 283.6 (l+Jcos2Trf r)e 2(5.66x10 ) 

p ^ in 

where 

T = (pAt) = p(20xl0 ® sec) 

f = 6.25 X 10® Hz 
m 

From Figure 9 it is clear that for dead times which are an 
appreciable fraction of the clock interval, there will be dis- 
tortion of the correlation functions. This does not appear to 
be a problem for At = Tm/8 except at the first delay valve and a 
general amplitude reduction but further study is needed. The dead 
time effect is seen to seriously affect the first delay; this 
will seriously affect the dual correlate approach if At = Tm/4. 

Variability Error 

In this section we derive formulas for the fractional rms 
error of a mean flow measurement with a dual correlate and sub- 
tract system in terms of the mean photon rate and the system 
response (integration) time. We assume that the accumulator 
sum is zero after summing for total time T. This implies that 
the clock frequency is in error by an amount required to cancel 
a random (finite time) variability error in the accumulator sum. 
These compensating errors are asstimed small so that a perturba- 
tion may be used. 
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Accumulator Sum 



An error in the accumulator sum is equivalent to an 
error e,p in the estimation of one quarter of the mean signal 
obtained by dividing by the slope of the discriminator function 
A(t) 




( 68 ) 


/N 

where A(t) is an abbreviation for <M > given earlier in equa- 

tion (56) and (57). The fractional error e in estimating the 

period is thus 
m 


e = 


V^'(V4) 4e^ 


T„/4 

m' 


T A’(T /4) 
m ' m' ' 


(69) 


This is the same form as equation (62) with the only difference 
being the type of error. We proceed by obtaining a simplified 
expression for A(Tjjj/4) by neglecting the Gaussian exponentials 
in equation (57). Evaluation of the derivatives gives 


A 


(V4) = 


NAt^R<X^> 

J 


4ir In 

f- V2 “ 

m 


(70) 


The quantity is the rms deviation of the accumulator value 
M after data collection time T. The evaluation of this 

pq 

quantity is discussed in Appendix C and the result for cases 
where the steady background light is larger than the signal is 
obtained.* For this case we have 


= «/ 2N 


(71) 


When these results are combined we have the fractional error 
in the velocity estimate given by 


*This is simpler to evaluate than the low background 
case and gives a bound on the time required to produce a given 
error in cases where the background is less. 
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c 


(72) 


o 

a/ TA t R<X .>tt/ it 

tJ 


where 


a 

<X^> 


R 

T 


background photoelectron rate >> signal rate 
burst 1/e half width 

mean square peak pedastal burst amplitude 
(photoelectron rate) 

rate of occurrence of bursts 

total time = NAt 


ADVANCED PHOTON PROCESSOR FOR 
TURBULENCE TIME STATISTICS 


Introduction 

In this section we propose and analyse a photon-processing 

scheme for estimation of the temporal autocorrelation of the 

time varying velocity fluctuations. Part of the basis for this 

is the fact that the frequency discriminator characteristic of 

Figure 5 applies not only to the long time average of the dual 

correlate and subtract sum 10 but also instantaneously in a 

pq 

conditional statistical sense. We may therefore tune At to the 
value which centers the long-term average at the zero of the 
frequency characteristic and then obtain short periodically 
occurring accumulations of Si whose conditional expected 

.r T. 

values follow the velocity deviation. Even when the photoelec- 
tron rate is small, the correlation of the velocity fluctuations 
may then be obtained by conventional digital correlation of 

the sequence of short time M ^’s. 

pq 

We define here AT as a period greater than At, the 
processor counting interval, and less than the significant 
times of change of the turbulent fluctuations. The quantity 
M is a sum over the interval t(n-l)AT < t < nAT] of m, 
defined in equation (55) as 
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(73) 


= “k<^k-p - ^k-q>’ 

i.e. , M =6 over the nth AT interval. We will consider 

pqn pq 

an autocorrelation estimate of the 0 sequence: 

pqn 

Rg(xAT) - ^pqn“pq(n+i) 

We will show that the expected value of this sequence contains 
the shape of the turbulence autocorrelation function under 
certain conditions. Further it appears that the magnitude of 
the turbulence intensity may be obtained by a normalization 
procedure which will be discussed. 

Spectral Analysis of Randomly Sampled Signals 

Continuous Funct ions . - In our previous development [23] 
we showed that correlations and frequency power spectra can 
be obtained from randomly timed discrete samples of a velocity 
component u(t) where the sampling function z(t) is a uniform 
Poisson impulse process. The same principles may be general- 
ized to the present more complicated data processing problem. 
First, let us assume that a continuous signal s(t) is available 
such that 


s(t ) 

= x(t)p(t) 

(75) 

p(t) 

= u(t)/u 

(76) 


where u(t) is a zero-mean time-varying velocity component 
deviation from the mean component U and the random process 
x(t) is a filtered Poisson shot noise process (see Appendix A) 
statistically independent of u(t): 
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x(t) = Ibf(t-T.) (77) 

j=oo J J 

Here f(t) is a low-pass impulse response function, the Tj's 
are random occurrence times which obey a stationary Poisson 
law, and the b.'s are equally-distributed statistically- 

J 

Independent random amplitude variables. The functions are 
illustrated in Figure 10. From the properties of Poisson 
processes we know that if f(t) is a positive function, then 
R (t) = <x(t)x(t + T)> is a positive function and, making use 
of the independence asstimption 

RpC^) = <p(t)p(t + T)> (78) 

= Rg(T)/R^(x) 

We observe that the zero value of R (x) is the normalized 

2 ^ 

mean-square turbulence intensity <p >. Figure 10 illustrates 

the fact that when the duration of the pedastal correlation is 

small compared with the duration of the velocity correlation, 

2 

the value <p > may be obtained approximately from Rp(x) where x 
is small but greater than the pedastal duration. If the turbu- 
lence intensity is obtained in some other manner, we may theo- 
retically obtain the shape of Rp(T) even without obtaining 
R (x), except in the vicinity of x = 0, since R (x) = constant 
elsewhere. 

We may extend the above reasoning to the situation where 
two sets of processes are available 

s^(t) = x^(t)p^(r^,t) (79) 

S2(t) = X2(t)p2(r2> t ) 

So long as x^ and X 2 are independent of p^^ and p^, then 
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Rs 12('^) = <Sj^(t)S2(t+T)> (80) 

= <Xj^(t)x2(t+T)xp^(rj^,t)p(r2,t+T)> 

and the cross correlation obtained by division 




(81) 


The frequency power spectra are Fourier transforms of the 
correlation functions. 

Sampled Data Estimates .- In the preceding, we assumed 
s(t) and x(t) were continuous and computed statistical corre- 
lation functions. We now suppose that real-valued discrete 
samples and of s(t) and x(t) are taken at uniform inter- 
vals AT which are small compared with the duration of the 
random sampling pulse f(t) and correlation time of p(t). We 
may form discrete estimates of R (t) and R (t) at x = iAT by 
computing a finite time average, for example 


Rg(iAT) 






n=l 


S S 
n n+x 


(82) 


where N^AT is the accumulation time T and the expected value 
is 


<R„ (iAT)> = R^(t = iAT) 

S S ' 


(83) 


If such estimates are formed for R and R_, then we may form 

X s 

an estimate for R (x) by division: 

p 


R (iAT) = R^(iAT)/R (iAT) 

p S X 


( 84 ) 


52 



This procedure does not insure that <Rp(iAT)> is an unbiased 
estimate but it becomes one in the limit of sufficiently long 
accumulation time when the numerator and denominator converge 
to their respective expected values. 

Conditional Expectation Again .- The problem is still 

more complicated. Because we do not have the classical optical 

signal X(t) available for direct observation, we cannot form 

x(t) and s(t), or even and S^. We will process the fast 

phot on- counting sequence nj^ to obtain discrete-valued sequences 

M _ and M whose conditional expected values given X(t) are 
xn sn c V 

proportional to X and . 

n n 

We thus have an estimate Rg(iAT) defined not by equation 
( 82 ) , but by 


The expected value of the estimate is 


<“sn“s(n+i)> = «“snl^(^)><“s(n+i)l^(^)>> 


= c <s .s . . .> 

3 J + i 


= c R (iAT) 
s 


2 

where C is the proportionality constant. This result depends 
on the fact that given X(t) the discrete random variables 
and formed from non overlapping portions of the 

sequence nj^ are conditionally independent. This follows from 
the fact that nj^ is Poisson counting process. 
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In equation (86) the notation <x|A> denotes "expected 
value of X given X" and we made use of the fact that the uncon- 
ditional expectation may be obtained in two steps: first, with 

respect to x given X ; then with respect to X . 

Application to Photon Processing of 
Turbulence Correlation 

It remains for us to identify measurable quantities with 
the properties attributed to and above and to describe 

a means of implementation. 

Conditional frequency discriminator .- We now reconsider 

/N 

the estimator M defined in equations (54) and (55). We 
pq 

restrict ourselves to the signal model given in equations (30) 
and (38) with X^ = 0 and with rarely overlapping bursts; 
i.e., a<<l/R. This is the "low density" shot noise case dis- 
cussed by Papoulis [lO] , page 574. For "low density" shot 
noise we may use the approximation that [10] if 


s(t) = Zh(t - T.) (87) 

3 

then 

s^(t) " 2h^(t - Tj) (88) 

with these constraints we obtain 

<nj^Hk-p I X(t)> = Ax^X(t)X(t - pAx) (89) 

oo 

= Ax^ X.^f^(t - x.){l + cos co.(t-x.) 

J=_oo ^ ^ J J 

+ cos w . ( t - X . - pAx ) + i cos(2o) .t - 2co . x . - pw . Ax ) 

J J ^ J J J J 

+ cos pwjAx} 


The Xi^ — 

the subtraction. 


0 assumption does not affect the result due to 
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In equation (89) the product of cosine terms has been expanded 

as sum and difference frequency terms, and the delay pAx has 

been neglected in the envelope function f(t). With mj^ defined 

by equation (55) we may determine <m, |X(t)> from two sets of 

terms like that in (89). We then determine <M > as 

pqn 


nL 


<M 


pqn 


|X(t)> = I <m, |X(t)> 
1 


(90) 


k=(n-l)L 


At 


nAT ( <» 

^n-l)4T 


- COS qtOjAx] + [Band pass term^j dt 


where 


L = AT/At (91) 

In the right hand side of this equation the sum is replaced by 
a time integral which it approximates. Since AT is long com- 
pared with the Doppler period, the bandpass terms average to a 
small value; however, since AT was assumed short compared with 
the burst envelope duration, the integral does not smooth the 
first expression in the right hand side of (90). The result 
is therefore 

00 

<M > “ y X f^(t - T . ) [cos pAxoi . - cos qAxco .] (92) 

^ j=-oo 3 3 3 3 

_ AtLT x^j^p(t)[cos pATw(t) - cos qAxaj(t)] 

By proper selection of a high speed digital clock period. 
At, the sequence approximates the sequence described 
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above for certain values of the integers p and q. For other 
choices of p and q the sequence ^pqjj approximates the sequence 
M^n* Thus, velocity correlations may be obtained by simul- 
taneously duplicating the computations of (73) with two differ 
ent sets of integers p and q, and then proceeding with soft- 
ware processing of the two sequences M and M . To show 

sn xn 

this last link in the procedure we examine the cosine differ- 
ence term in (92). For the M^^sequence we select p and q and 
adjust the variable At to the conditions given previously for 
the mean frequency discriminator in equations (58) and (59). 


(P - (93) 

PO)„At = 37r/2 
m 

qw^Ar = tt/ 2 
m 


Example selections for p, q and At are 3,1,T /4; 6,2,T^/8; 
12, 4, T 716; where = 2 tt/w_ is inverse of the mean signal 
frequency. Expansion of the bracketed cosine term in (92) 
gives 


[cos pAto) . - cos qATO) 


ttAo) . ottAo) . ' 

m m -J 


(94) 


where we have used the substitutions 


Au . = 0) . - 0) 

J J m 


(95) 


cos (Wj^ + Auj) = sin(^^j^ — i) 



37rAu. 
''m 


& 


l2Zr <“m * 

m 


= -si 


irAo) . 
m 
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Equation (94) was shown plotted in Figure 5 since it is also 
the approximate long-time average frequency transfer charac- 
teristic for the mean flow measurement. The trigonometric 
substitutions along with the small angle approximation sin 0-0 
show that for small turbulence levels, equation (92) becomes 

<Msn|X(t)> = ATAtX^j^P (t)Trp(t) (96) 

Where p(t) » non-zero portions of ^j^p(t)* 

Normalizing sequence .- Under the same selection of At as 
in (93) but with 


pu^Ax = 2 tt 
m 

(97) 

> 

II 



a different effect is obtained. For example, if M is obtained 

sn 

with p=3, q=l and we let p = 4, q*2 in computing M , the 

All 

result corresponding to (94) is 


[cos pAxo)^ - cos qAxw .] 



(98) 


The small-angle approximation for the cosine function is unity, 
so the result corresponding to (96) is 

<Mxn|X(t)> = ATAxX^^p(t) (99) 

When we note that Ato./w^ = u./U and that the ATAt term cancels 

3 m J 

by division, we see from equations (96), (99) and (86) that 

the autocorrelation of M„ divided by the autocorrelation of M „ 

sn •' _2 xn 

produces approximately ir <u(t )u(t + x )>/ U . Except for a fac- 
2 2 

tor of C = IT , the result is normalized in such a manner that 
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fractional mean-square turbulence intensity is directly obtain- 
able from the first lag value of the final autocorrelation 
estimate or a projection back to the zero delay value. (The 
analysis we have presented is not valid for the zero delay 
value without separate consideration. It is possible that cer- 
tain effects will cancel and cause this to be a valid point 
also. ) 

Practical considerations .- The approach outlined above 
for normalization would require as much hardware to measure the 
sequence as the sequence. In addition, the small angle 
approximation which was applied to (98) is not valid over a 
very large range of velocity deviation. We have conceived two 
other less expensive approaches. The first of these consists 
of eliminating the second M channel completely and relying on 

aXX 

the fact that the shape of the correlation Rp(T) is approximated 
by correlating except in the vicinity of a burst duration 
from the delay origin. If it were desirable to normalize the 
function at the rest of the delay locations, we could evaluate 
the required division constant as 

AT^At^^ = <Rg(nAT)>, nAT >> a (100) 

2 

where we have already discussed a possible estimate for <A (t)> 
in equation (46). Either this estimator or the approach 
described above could be computed with a single channel elec- 
tronic system by electronically switching the delay variables 
p and q and storing the sequences in different portions 

of a computer memory before computing the second level corre- 
lations. Time sharing like this would not be as efficient in 
possible cancellation of some of the statistical variability 
error, but significant cost savings would result. 
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DISCUSSION 


System Design 

We have provided a specification for an advanced computer 
controlled photon counting processor in Appendix D. The design 
allows the system to be operated as a dual-correlate and sub- 
tract processor or sequentially as a correlator by time sharing 
the multiplier. The system provides synchronous 3-bit X 3-bit 
operations at up to 100 MHz. The design uses "slow" emitter 
coupled logic (ECL) which is optimum for the speed range speci- 
fied. The system may be operated as an advanced processor by 
dumping the accumulator values to Computer memory at rates 
limited only by the computer. The design allows two identical 
units to be used together for either velocity cross-correlation 
measurements or for simultaneous second-channel normalization. 

An analog feedback loop is provided for automatic mean velocity 
acquisition. In addition, the system can be used to measure any 
type of multiple-interval photon statistics by selecting the 
correlator clock period (continuously variable both manually and 
electronically) and using the computer to sample the values of 
stored in the delay line. 

Sensitivity Comparison 

There are four primary sources of variability error. 

These are the random turbulent flow itself, the random occur- 
rence times of the scattering particles, the random scattering 
cross sections, and the random time of photon events. Photon 
correlation methods are linear in the sense that the effects of 
two simultaneously occurring scatterers are added. This is 
beneficial in that it avoids the non-linear zero-crossing cap- 
ture effect of classical FM systems and thus the error problems 
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of multiple scatters; it is not beneficial in the sense that 
random amplitude effects will contribute to variability error, 
but these distinctions have not been carefully analyzed. The 
error versus data collection time trade-offs which result due 
to the random occurrence times of scatterers and due to mini- 
mum limits due to the turbulence itself have been previously 
analyzed for burst counters. That result* may be expressed as 


T 


<P^> 

<e^< 



( 101 ) 


where B is the equivalent power bandwidth of the turbulence, 

2 

R’ is the mean rate of accepted signal bursts, and <p > and 
2 

<e > are the normalized mean-square turbulence intensity and 
mean-square estimate error, respectively: 


<p2> = 


( 102 ) 


<^2> = <ll - U)^_> 


this result indicates, for example, that for a turbulence equiva- 
lent power bandwidth of 2KHz, 10% turbulence intensity, and 0.3% 
desired rms error, the data collection time would be at least 
0.28 seconds even if the continuous signal U(t) were available 
for processing, and would be considerably longer if R' were 
less than 4000/sec. 

In the following we examine the implications of equation 
(72) for transonic wind tunnel measurements. In order to include 
the mean value of the signal bursts along with the assumed 


page 10, equation (2.14), reference [23]. 
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background light in equation (72) we use equation (40) for 
<X(t)> and (72) becomes 


e ~ 


2(X- +X ) 
' b s'^ 


u/TAt R<X •> Ti /it 
3 


(103) 


where the mean rate due to signal bursts is 


X = R<X . > /ir a 
s J 


(104) 


which is less than the mean peak rate when R is less than 

the inverse of the effective burst width, /ir a. Now by defin- 

2 

ing the variance a, of the pedastal peak rates X . as 

A J 

o ^ = <xj> - <X.>^ 

X J J 


we may rearrange (103) as 


e = 



(105) 


Although the parametric behavior of equation (105) is 

intuitively acceptable in other ways, the presence of the term 
2 2 

[1+0, /<X.> ] in the denominator seems a little strange. 

^ J 

Increased variability of the scattering cross sections of the 
particles would not intuitively decrease the mean-flow varia- 
bility error. It is possible that this indicated behavior is 
a consequence of ignoring the variability of the classical power 
X(t) in equation (71) and has no physical reality. However, 
there does not appear to be anything wrong with the derivation 
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for the limiting case of high background light so the strange 

result may be true. In any event, a conservative bound is 

obtained by removing the bracketed expression in the denomina- 

2 

tor to obtain an expression not requiring knowledge of : 



(106) 


If we evaluate equation (106) with the following assumed typical 
conditions for a transonic wind tunnel measurement we obtain 
a required data collection time of 0.5 seconds for a 1% rms 
error : 


\ = 10*^ (107) 

<X.> = lo”^ 

u 

At = 10"® 

In this example, the mean signal photoelectron rate is ten 
times less than the mean background photo-electron rate and 
is equal to the average peak envelope rate. With this much 
background light, the assumption of constant X(t) should be 
valid with respect to the photon-fluctuation induced varia- 

7 

bility. The selection of mean peak rate at 10 means that 

g 

occurrences of photoelectron count rates greater than 10 /sec 

(the limit of current hardware state of the art) will be rare 

and the effects of nonlinearity negligible). The selected 

ratio of <X.>/X = 10 implies that the measurement volume is 

j s 

only assumed to contain a scatterer 1/10 of the time on the 
average. 
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Even though dramatic improvements over the previous exam- 
ple situation may result from reduced background light, the 

result still indicates that practical measurements may be 

7 

obtained with a peak signal photoelectron rate of 10 sec. In 
order to compare this with the performance of a burst-counter 
system, we must assume values for p, q and f^^^, the Doppler fre- 
quency. With p = 3, and q = 1 the Doppler frequency f^^ is 
1/4 At = 25 MHz. This would result from U = 30^.8 m/sec with 
an optical sensitivity of 82 kHz /m/sec. This peak photoelec- 
tron rate assumed is thus 0.4 photoelectrons/cycle in the pres- 
ence of 0.4 background photoelectrons/cycle. For comparison we 
note that standard optical noise formulas, given X(t) = 

Xj(l + cos w^t) + X^, would result in peak SNR of 


SNR = 




4B(Xj + 


(108) 


7 9 

In our example, B = 25 MHz and X, = 10 ; if we choose X . = 10 , 

® J 

a factor of 100 greater than in our example, then the SNR is 
10 at the peak of the signal burst and 1.35 at the 1/e signal 
envelope points. Since this example represents marginal or 
inadequate SNR for burst counter operation we deduce that, even 
with 100 times more scattered power, only the larger-than-aver- 
age scatterers would contribute. 

Under conditions of less background light, the burst- 
counter analysis would not be improved; however, the photon 
counting system results are expected to improve considerably. 
Thus we conclude that mean-flow measurements with from 100-1000 
times less optical power are feasible with the photon counting 
system. 
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Future Work 

During this contract effort we have developed simulation 
techniques which can be extended to be applicable to all types 
of LV signal processors for any level signals. The program in 
Appendix B only represents the beginning of how these techniques 
can be exploited, and we did not have time under the present 
contract to use that program except for check-out waveform simu- 
lations and the dead-time example. Similarly, the higher order 
moment equations developed in Appendix A have not been yet used 
to extend the variability error analysis to include the low 
background case. 

A system such as that specified in Appendix D should be 
constructed and tested. The results in this report indicate 
that such a system would a low LV measurements to be performed 
which are not now feasible. In addition, the system would be 
a valuable research tool for many other fields of research 
which require high-speed digital correlation or measurement of 
photon interval statistics. 
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CONCLUSIONS 


The most significant new results which this report provides 
are itemized below: 


• A general non-st at ionary Poisson process model for 
dual sc,atter laser velocimeter (LV) signals and noise 
valid from high level signals down through photon 
resolved signals. 

• Computer simulation algorithms valid over the entire 
range of signal levels, which may be used to evaluate 
any new type of LV signal processor. 

• A description and statistical analysis of both con- 
ventional photon correlation and Dual Correlate and 
Subtract frequency discriminator technique for mean, 
turbulence intensity and turbulence correlation 
estimations from photon resolved signals. 

• A system design for an advanced photon-counting 
processor which implements both conventional photon 
correlation (sequentially) and the Dual Correlate 
concepts with time resolution to 10 nsec. 
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APPENDIX A 

FILTERED INHOMOGENEOUS POISSON PROCESSES 

In this appendix we have used the material from Snyder's 
book [9] to derive the higher order moments of a filtered 
inhomogeneous Poisson process up through order 4. We begin 
with introductory material similar to that given by Papoulis 
[ 10 ]. 

Inhomogeneous Impulse Processes 

The input to a random linear system is an inhomogeneous 
Poisson impulse process z(t) given by 

oo 

z(t ) = I 6(t - T^) (Al) 

— OO 

where is the set of random occurrence times, X(t) is the 

instantaneous statistical mean value of z(t), (and also the 
mean rate of occurrence of the t^'s), and 6(t) is the dirac 
delta function. The random variables are independent of 
each other statistically and obey the inhomogeneous counting 
law, i.e., the probability of n=k occurrences in the interval 
(t^,t 2 > is 

_y 

P{n(t^,t2) = k} = <A2) 

where 

^2 

y = / A(t)dt (A3) 

h 

The quantity y is also the mean and variance of the random var- 
iable n(tj^,tg). 
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The Response of a Random Linear System-Campbell' s Theorem 

The output of the random linear system s(t) is the super- 
position of the response h(t - Y^) to each input impulse: 

oo 

s(t) = I h(t - Y^) (A4) 

00 

where is a set of identically distributed independent 

vector random variables. The random variable Y affects the 
shape and amplitude of the response function h(t,Y). In the 
case of the PMT signal it may take the form of a single scalar 
amplitude variable. In the case of the classical optical sig- 
nal from turbulent flow both a random amplitude parameter and 
one or more random shape parameters due to velocity magnitude, 
direction, and probe volume translational entrance location may 
be required. The theory should be applicable so long as the 
set of multidimensional random variables Y^ is independent of 
the set of occurrence times {t^}. The generalized Campbell's 
theorem results for the instantaneous statistical mean, variance, 
and auto-covariance of s(t) are given below, they apply regard- 
less of whether individual pulses are resolved or not. 

OO 

<s(t)> = / X(t ) <h(t - T , Y)>dx (A5) 

->oo 

oo 

o^(t) =<s^(t)>- <s(t)>^ = / X(t ) <h^(t - X , Y)>dx (A6) 

OO 

cov [s(tj^)s(t 2 >] = <s(t^)s(t 2 >> - <s(t^)xs(t 2 )> (A7) 

OO 

= / X(x)<h(tj^ - x,Y)h(t 2 - x,Y)>dx 

OO 

where <> denotes expectation with respect to Y inside the 
integral signs. For a causal signal such as that from the PMT, 
where h(t) is zero for t >0, then the upper limits of 
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integration may be replaced by t or the minimum of t^ and tg 
in equation (A7). For the transient case where the impulse 
signal z(t) is applied at t = 0, the lower limits of integra- 
tion may be replaced by 0. 

Higher Order Moments 

Summary of statistics .- Given a filtered, inhomogeneous, 
compound Poisson process: 


oo 

s(t) = I h(t,T.;Y.) 

— OO 


(A8) 


where Xj^are random occurrence times which occur with intensity 
X(t), and where are vector random variables which are sta- 
tistically independent and identically distributed, we obtain 
the following result. The cumulants are: 

OO 

Yi(ti) = = / X(t) <h(t^,x;Y) > dx (A9) 

OO 

oo 

"*'l2^^1’^2^ = ^ / A(x) <h(tj^,x;Y)h(t 2 ,T;Y)>dx 

OO 

oo 

^123^^1'*2’*3^ ^ A(x) <h(tj^,x;Y)h(t 2 >x;Y)h(tg,x;Y)>dx 

oo 

oo 

'^1234^*1’^2’^3’*4^ = / X(x)<h(tj^,x;Y)h(t2,T;Y) 

00 

• h(t 3 ,x;Y)h(t^,x;Y)> dx 


The formulas which relate the cumulants to the moments 
are as follows: 


1 . 

2 . 


<s(t^)> = = Y;L 


<s(t^)s(t2>> = 

where y ^2 ~ covariance = 


(AlO) 
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3. <S(t^)s(t2)s(t3>> = Yi23 + ^1^23 "^2^13 

■*■ ’^3^12 ’^1^2’^3 

4. <s(tj)s(t2)s(t3)s(t4>> = Y]^234 ^l’^2^3’^4 

+ ’^ 1^234 ■‘' ’^ 2’^134 ■‘' ^ 3‘^124 ’^ 4^123 

■*■ ^12^34 ■*' ^13*^24 ■'' ^14^23 
■‘‘ ^ 14 ^^ 2^3 ^ 13 ^ 2^4 ^ 12 ’^ 3^^4 

•" ^24^1^3 "■ ^23^l’^4 ^ ^34^1^2 

The derivation of the above formulas follows. 

Derivation . - The derivation of the preceding formulas is 
straightforward but tedious if we are given the joint charac- 
teristic function 

. 03 : s 

^(( d ) - <e'^ > (All) 

where 


and 


Cl) (cij^^j ^2’ ^3^ ^^3 


(A12) 


s = [s(t^),s(t 2 ),s(t 3 ),s(t^)] (A13) 

From Papoulis we know that 
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(A14) 


<s(t^)s(t2S(tg)s(t^)> 


1 3"^(j)(to) 

(J)^ StOg BWg 3(1)^ 


w=0 


and similarly for other moments; i.e., we may obtain the moments 
by determining the appropriate partial derivatives of the Joint 
characteristic function at co = 0 which is in turn a vector 
valued Fourier transform of the joint probability function for 
the random variables s(t^), s(t 2 ), etc. 

The theory of filtered Poisson processes provides us with 
the joint characteristic function of the second kind, ip, where 


ip(w) = in (p (to) 

(A15) 


(A16) 


We may therefore either write out (p (w) and evaluate the 
derivatives directly or use the above equations to first 
express the partial derivatives of <j) (to) in terms of the par- 
tial derivatives of t(j(to). We have taken the latter approach 
using the chain rule. As an example 


3 Wj ^ 3(02 ^ ' 3 ( 0 j ^ 3 t 02 ^^1 ^^^2 ^ 


(A17) 


The rest of the derivatives are omitted here since they get 
progressively more lengthy. 

We now need only to evaluate the products of partial 
derivatives we have obtained at to = 0. From the material given 
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1> dT 


(A18) 


\ 

f 


by Snyder we have 

- 

— oo 

where 

3(w) = w^h(t^,T;Y) + UghCtg.x;?) + . . .+ o3^h(t^,T;Y) (A19) 

From the form of the expression for i|^(aj) it is clear that 

I “ / X(x)<e^^ - 1> dx = 0 (A20) 

I — _oo 

o )=0 

therefore 


4>(w) I =1 

( i )=0 


(A21) 


It is not difficult to see that the form of the partial 
derivatives is 




oo 

= / X(x) <h(t. ,x;Y)h(t 5 ,x;Y). .e'^^^“^>dt 

. I w=0 -» 12 


( i )=0 

(A22) 


i.e., the partials of 4^(o)) at u = 0 are equal to the same order 
cumulant except for the factor. With this formula we can 
now go directly from the expression of the partials of <p in 
terms of the partials of ip to the desired higher order moments. 
The result is that given in the summary. If moments higher 


♦ 

In Snyder's book, Random Point Processes , equation (4.15) 
the lower limit of integration is to correspond to the begin- 
ning of the process. The upper limit is the minimum of the 
times t.. This assumes that h(t.,x;T) is a causal function. 
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than the 4th are needed we would simply apply the chain rule 
to determine the higher derivative of <f> in terms of those of 
Ip. 

Moments of a Gaussian Random Process .- We have not yet 
expanded all of the moments of a non-zero mean Gaussian 
process for comparison. We observe, however, that the factori- 
zation property of zero-mean Gaussian processes does not hold 
for filtered Poisson process. Otherwise there are great simi- 
larities except for the cumulants as 


Gaussian 
<xi> = n-^ 

<XiX2> = ’^l’^2 

for (n^ = 0) 

<XiX2Xs> = 0 

<XiX2XgX4> = + M13U24 

* ‘'14'‘23 


Poisson (A23) 

<Xi> = rij^ 

<XjX2> = '^l'^2 

<XiX2X3> = Yi23 
<XiX2X3X4> = y^2^34 ^13^^24 

^14^23 *’^1234 


(n f 0) 

4 4 2 2 4 

<x^> = 3a + 6a n + n 


<x^^> 


3a^^ + 6n^a^ + 


4 

n 


+ 4n / Xh^(t)dt 
+ JXh^(t)dt 
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APPENDIX B 


DIGITAL SIMULATION OF LOW LEVEL DUAL SCATTER 
LV SIGNALS AND IDEALIZED PHOTON PROCESSORS 

This appendix provides theory and an example FORTRAN 
program for digital LV signal simulation. Background theory 
and several more complicated Poisson impulse simulation algo- 
rithms are discussed by R. H. Forrester, Jr. in a Masters Thesis 
[25] performed under Don Snyder. The simple approach taken 
here is to discretize the possible occurrence times of photon 
and classical burst signal events to uniformly spaced inter- 
vals of length At, where At is less than any significant sys- 
tem integration time. This discretization of time imposes 
itself upon all temporal system parameters, such as processor 
clock interval At and dead time t^, which would be continuous 
variables in the real situation. 

Theory 

Realization of Inhomogeneous Poisson Impulse processes . - A 
sample function of an inhomogeneous Poisson impulse process is 
specified by a set of event occurrence times {t^} as described 
in Appendix A. The key procedure required for simulation of 
LV signals is therefore the generation of a set {t^} given a 
specified rate intensity function X(t). This may be done by 
first generating a homogeneous (stationary) process with unit 
intensity X(t) = 1 and then mapping the realization times 
through solution of an integral equation. Figure B1 illus- 
trates the required mapping in a manner which helps provide 
intuitive grasp for what follows. 

Let us define the set of interarrival intervals {w^} 
between event occurrence times as 
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unit intensit 



w. 

1 


(Bl) 


For a homogeneous Poisson impulse process with constant rate 
intensity X(t) = X, it is necessary and sufficient that {w^^} 
be a set of statistically independent, identically distributed 
exponential random variables with common probability density 
P^(w), given by 

(B2) 

We may therefore generate the set ' of occurrence times for 

a realization of a constant intensity process as 


T . 
1 


W, 


= X 


k=l 


i-1 


+ w, 


(B3) 


Commonly available subroutines generate statistically 
independent realizations of a random variable x uniformly ' 
distributed on the interval (0,1). The exponential density 
function is monotonic with an inverse function which is commonly 
available (namely the natural logarithm function). It is 
therefore straightforward to determine a transformation which 
maps the realizations of the unit uniformly distributed varia- . 
ble to the desired exponentially distributed variable: 


”1 = 1 

In summary, generation of unit uniform random variables x; 
which are then transformed by (B4) and used in (B3) produce a 


♦ 

There is, however, a difference in the quality of these 
subroutines; a survey paper by [26] may be consulted if the 
validity or efficiency of a subroutine is in question. 
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realization of a homogeneous impulse process with intensity 
X. For our purposes, we set X = 1 and proceed to the mapping 
step illustrated in Figure Bl. 

In his recent book* Snyder provides guidance, but leaves 
it as a homework exercise to deduce the proof for a method of 
rescaling the interarrival times of a unit intensity process 
to obtain a realization of a specified inhomogeneous process. 
The results are as follows: Let be the set of occurrence 

times of a unit intensity process as illustrated in Figure Bl. 
Let A(t) be the integral of the specified intensity function 
X(t): 

t 

A(t) = / X(a)da (B5) 

o 

This function is continuous and monotonically Increasing and 
therefore has an inverse function A”^ such that 


- A(t^) 

(B6) 


(B7) 


The set of times {t^} generated by applying equation (B7) to 
the set is the required realization of the inhomogeneous 

process. 

Algorithm for simulation of event times .- The subroutines 
which generate uniform random variables produce real numbers 
If the above equations are applied exactly, one would be forced 
to digitally solve an integral equation, thus producing another 
set of real numbers {t^}. One could carry the full resolution 


♦ 


See Snyder [9] page 62. 


76 


of the computer to the bitter end of the simulation. However, 
at some point the use of digital filtering techniques to 
simulate analog filter characteristics of real photomultiplier 
tubes (PMT) and other electronic devices would become appro- 
priate. At such a point, the simulated signal would have to 
be interpolated and respecified on uniformly spaced increments 
with separation At. We may greatly simplify the required 
algorithms by rounding the set of occurrence times to the near- 
est At interval and point. By incorporating this step directly 
into the solution of the inverse function we avoid 

the problems of solving the integral equation exactly and sim- 
plify that step as well. The entire procedure is thus simpli- 
fied to the following. 

1. Select a At small enough to provide adequate 
accuracy for uniform sampling of X(t) and cal- 
culating its integral by the trapazoidal rule 
integration method. 

2. Beginning at ” ^o ” compute realizations 
of as discussed above. 

3. Calculate the trapazoidal rule approximation 
of A(kAt) for each integer value. ® — ^ — ^max' 

4. Use conditional statements to test the latest 
real value of against the integral A(kAt) 

as it is generated iteratively to determine a 
histogram TC of the discretized occurrence 
times at kAt . It is possible, in the simula- 
tion, for more than one value of x. to be 
mapped to the discretized time kAt. 

There are many ways that the fourth step could be implemented. 
We have elected the following: If the value of x^ lies in the 

range 

A(kAt) + A[(k-l)At] < < A[(k+l)At] + A(kAt) 

2 i 2 

where A (kAt) is the trapazoidal approximation of the integral. 
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then an occurrence is added to the histogram value at k. 


Inclusion of random amplitude effects .- When the point 
event represents a photoelectron pulse or a classical LV sig- 
nal burst, we may wish to assign the event an amplitude mark 
(impulse weight). This may be done by separately generating 
realizations of additional random variables according to 
desired statistics and accximulating one each of these at TC 
for each occurrence time. Note that it is not adequate to 
simply multiply each value TC by a random variable, since in 
some realizations more than one event contributes to the same 
value of k. 

The statistics of the pulse height distribution of a poor 
(PMT) may be nearly Rayleigh, while the Gaussian density with 
15-25% relative standard deviation may bo adequate to model a 
good PMT. Very little documentation exists concerning the 
probability density of the classical signal bursts. This topic 
is discussed and some data is presented in our recent AEDC 
report [1] . Of the simplest densities an exponential density 
or a Rayleigh density would be used to simulate amplitudes from 
an unseeded flow. As we have shown it does not follow that 
the amplitude probability density agrees in any recognizable way 
with the particle size distribution; even monosized particles 
may produce a very strange amplitude probability density [1] . 

When it is desirable to use Rayleigh or Gaussian random 
variables the following procedure is recommended. Generation 
of realizations of Rayleigh or Gaussian random variables may be 
obtained by first generating uniformly distributed values on 
the interval (0,1). Let x^^ and Xg be two such independent 
realizations. Then we obtain 

R = (-2a^ In (B9) 

0 = 2-nx^ (BIO) 
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where R is a realization of a Rayleigh random variable with 
parameter a, mean = , and probability density Pjj(R) given as 


Pr(R) 




(BID 


and 6 is a realization of a uniformly distributed random vari- . 
able on the interval (0,2 tt). Multiplying (B9) by / 2/ it produces • 
Rayleigh variables with mean =1. If Gaussian random variables 
are desired, the process is continued from (B9) without the 
/2/ir factor by converting to rectangular co-ordinates: 

(B12) 
(B13) 


X = R cos 6 
Y = R sin 6 


When this is done, X and Y are two independent realizations of 

2 

a Gaussian random variable with zero mean and variance a , i. e., 

P„(X) = (B14) 

and the same form for Py(Y). The above procedure efficiently 
produces exactly Gaussian random realizations as opposed to a 
program such as GAUSS which sums 12 independent uniform random 
variables to obtain approximately Gaussian variables by the 
central limit theorem (Forrester [25]). 


Example Simulation Program 

We have included at the end of this appendix a copy of 
the printout of a FORTRAN IV program which is illustrated in 
flow form in Figure B2. The occurrence times of the classical 
signal bursts are generated as a homogeneous Poisson process. 
The amplitudes of the bursts may be either generated randomly 
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Figure B2. 


Flowchart Showing Namelist 
Names and Functions. 


Variable 
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with exponential, Gaussian, or Rayleigh density with a speci- 
fied mean, or they may be set equal to a constant. The pro- 
gram computes the burst waveforms at each At and sums all that 
are present to form a classical signal. This signal is inte- 
grated in a trapazoidal fashion and compared with the occur- 
rence times of a unit intensity homogeneous Poisson impulse 
process simulated as described previously. If random amplitudes 
have been assigned to the photoelectron pulses, these are gen- 
erated and added to a histogram; otherwise I's are added to the 
histogram. 

The idealized photon processor portion of the program sums 
the values of the histogram (with no random PMT pulse height 
effects) over an interval At which is some selected integral 
number ( ITAU) of At units in length. The sum is the photon 
count sequence referred to in the text as The ideal 

processor then computes either the photon correlation for delay 
values 0-IP (ONE = false) or computes and sums the dual corre- 
late and subtract terms (ONE = true). 
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: 83818 F 33 ISI FtBft 833 HBE 

= logical rHICH uses ideal processor RHEn set true 
= LIMIT OF THE COUNTER IN IHE IDEAL PROCESSOR 

' |pl8{fiB 1^‘^AHl/irA^^Y'^ibVpuT file for PLO.S 
= THE 1/E Squared half riuih (seconds) of bursi 
= determines RhERE burst la TRUNCATED TO ZERO 

= .^i58JekTfir?N=o^ra^ BUHSTSCE**(-2T**2/0**2)* 

= the integer number of OT a COMPRISING ONE PROCESSOR 
OTAO 

= THE REMAImOER of AN INTEGlR DIVISION OF THE NuMBER 
OF DT S BY ITAO 

s SECOND CORRELATION DELAY INDEX AS USED IN 

s cS^^^aI ^0^^3eL^Y {n§^^ maximum or FIRST VALUt IN 
DUAL correlate 

s SUM OF PHJTUELECTKUNS COUNTED DURING TOTAL MEASUREMENT 
= accumulator sum 

= logical RHICH IS aE I TRUE WHEN A CLASSICAL BURST IS 
present 

S logical which is SEI TRUE WHEN A CONSTANT BACKGROUND 
IS PRESENT 
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1600* UTIML 
1610* IPOJN 


1 630*C 
1640*C 
l650*C 
1660*C 
1670*C 
1660*1 
1690*C 
1700 * 
1710*- 
1720*- 

175>0*C 

um 

1780* 
17Bb* 
1790* 
>* 
i610* 
16?0* 
1630*C 

{elo*c 

1660 * 

1 670*L 
1880*C 
1690* 

1 900* 

1910*C 

I920*t 

1930*1 

19/10* 

1950* 

1960* 

1970* 

1 960* 

1 990* 
2000*C 
2010*0 
2020* 
2030* 
20/t0* 
20/15*0 
20/l6*C 
20/!7*0 
2050* 
2060* 
2070* 
2060*0 
2065*0 
2066*0 
2090* 
2100* 
2110*C 
2120* 
2130* 
2140* 
2145* 
2150* 
2160* 
2170* 
2160* 
2165* 
2190* 

2220 * 

2230* 

2240* 

2250 * 

2260 *- 


- CONS 
Dt-AD 

- 

" ot 

NAMFl I 
D.CPN 


ST 


w, KMAX 
IP^PAXC ‘ 


s THF constant background LAMBDA TO BE ADDEO IN THE TRUE 
CUNDITION or const 

« THE integer number UF DT S OF THE DEAU TIME INTERVAL 
IF NOT USED then IT DEFAULTS TO ZERO. , . 

= LOGICAL where TRUE DUAL CORRELATE/ FALSE„ I MPL I ES„,, 

IP IS the max delay in PHOTON CORRELATION, 

« pRoriiCT OF electronic charge and charge gain g for PMT 

/INP/If.l T/ I nAUX/ THE OR Y/ plots, PHI/ PLO/ WAVE, CONST /CONS/ ONE/ 
MAX/Tb, Ar AC/ HIGH/ LOW/ I R 1 / 61 / 1 R2/ Ci2, 1 Q, GE, AME, LO/ 1 0 A, 

/DEAD/ iTAll/TTUT/MNTOT/TUTlME/NTOT/ptPME/MTOT/DT 

1Ki PPriRLAM UItH a KlAMriT.<i1 DFAO 


- HFAD IN PRDGFAM VaLULS WITH A NAMELIST READ 

m 


[ ^C = f/j/D/CO^'*KMAX 
DT=D*cnN7KMAX 
1DTT0T=TTDT/DT 
UTIMEsRANDC 1/ JV/G2)*F AC 


BURSTSd Of 
PEAK lambda 
BOUNDS 


burst LOOPS) 

generation 

( HIGH/ LOW )*AFAC 
TO 15 


. BURSTS 
bEl WEEN 


= 0 

P0JNT=1 
Mb 4 

pEPlORM A SIMUIATION FOR N$B SIGNAL 
999 CONTINUE 

ftt-T A RANDOM VALUE FOR THE PEDaSTAL 
15 AMPbRaNU( 1R2, JV/GI )*AFAC 

- IF THIS value IS NOT IN PEAK LAMRDA 

- THEN GO BACK AFD TrY AGAIN 
1F( AMP.GT , ( AF AC*H]GH) ,0R, AMP.LT , ( aFAC*LOW))GO 
JENDsUTIME-OTIFF 

- DETEKMINE^IF there IS SPACE IN BETWEEN THE SIGNAL 

- ****THIS BRANCI IS USED FOR SPACING OR FILLER IN 

- sicnal Bursts**** 

lFr.lENI),LE.O)GL TO 3& 

IFC .NUT .PLOTS, CP .CONST )G 0 To 99 
L>045J=1/<JFNO 

45 CALL pLUTSpd / J/KMAX/0T1 mE/0,C0N/ JOAUX) 

99 OTIMEbOT IME+JFNO 
/lEROSc.TKUt , 

- THE FOLLOWING FROCfDURE (Up TO STATEMENT 46) IS 

- REPUIPEU for CASES THAT HAVE A CONSTANT BACKGROUND 
IF( .NUT .CONST )cr TD 46 

1 ESTbWAVE 
WAVFb.f ALSE. 

- here we evaluate the space in between signal BURSTS AT 

- INTERVALS OF A 1000 DTS AT A IIME/UNTIL THE SEpLRATlON IS 

- <1000. THEN WE JUMP TO 47 AnD EVALUATE NORMAL SePERAIION. 

50 lOLApBjF ND-1 ( 00 
IrC T0LAP)47/ 47//<9 
49 7 EH()Ss,1RUE. 

- CALL t)N the SK-S ROUTINE UnOER THE ZEROS CONDITION 

- because max afpay size is 1000 WE evaluate for MAXIMUMS 

- OF 1000 DTS AT A TIME, 

call S1GS( AMP/ ICOO/ 10/ IRl, I r 2/D,CDN/ AME/F0/WAVE,C0NST/C0NS) 
-^CALL -6S‘'r^E ’ ICf AL PROCESSOR 

1F( THEORY) Call IDEALC 1000/ IPOINT/ it AU/ IOA/ IP/NT0T,MTUT/ IR/DEAD) 

JENDbIOLAP 

Uo 51 J=l/1000 

IF(PLOTS) CALL F I OT ( M/ J/ KM AX , OT 1 ME/ D/ c ON/ 1 0 AUX ) 

R( J>=0 
TCR( J)=0 
S(J)=0 

51 TC(J)=0 

IF(DONE) go to 996 
Gn to 50 , 

47 ZER0S=.TRUE. , , 

- CALL ON THE SIGS ROUTINE UNDER THE ZEROS CONDITION 

call SIGSCAMP/ JFND/Iti/IRI/ Ir2/U,C0N/AME/E0/WAVE/C0NS1 / CONS) 
WAVE=TEST ^ 

ZER0S=. FALSE, 

46 IF (THEnRY)CALL I OE AL( JEND/ I PO I NT/ I T AO/ 1 0 A/ I P/ NTUT/ 

MTOT, IR,UEAD) 


S3 


?270* 

2260 * 

2290 * 

2295* 

2300* 

2310* 

2320* 

23^tO*C« 

2350* 

2355* 

11 ^ 8 : 

2380*C* 

2390*C* 

2 « 00 * 

2A10* 

2420* 

2430*- 

2431* 

2432* 

2433* 

2434* 

2435* 

2440* 

2450* 

?Hf 8 : 

2460* 

2490* 

2500* 

2520* 

2560* 

2570*C* 

2560*C* 

2590* 

2600 *- 

2610 * 

2630* 

2640* 

2650* 

2655* 

2660 * 

2670* 

2b7b* 

2660 * 

2690* 

2700* 

2710* 

27 ? 0 * 

2730* 

2740* 

2750* 

2760* 

2770* 

2760* 

2790* 

2600 * 

2610 * 

2620* 

2630* 

28<'0* 

2650* 

2655* 

2660* 

nm 

2690* ’ 
2695*- 
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ZERHSs, FALSE. 

IF( .NOT.CONSncr TO 35 
On 46 Jsl^JENO 

IF(PL0TS)CALL FI 0T(M* J,KMAX,OTIME^O*CON#IOAUX) 

R(J5=0 

TCR(J)s0 

S(J)=0 

^CALL^ 0 N~?HF sir-s routine to compute a signal burst 

35 CALL SIGS(AFP.KMaX,I0,IR1/IR?/D/C0N^AME#F0,WaVe#C0NST#C0NS) 
IF(PONE) GO TO. 996 ^ ^ , 

lF((KMAX*oHME)CL?iDTfME)Gi^^O 25 

**** this brafch is used for The overlap condition 
OF signal bursts**** 

INTFR=DTIML-0TU'E 

lNSECT=,TRUEf 

If (IHEORYICALL I DE aL ( I NtER. I PO l NT* I T AU* I OA/ I P# NT OT , HTOT* 
IR.OEAO) 

DO 41 jsl.iNTEF 

IK PLOTS KAIL F L OT ( M* J* KmAX, OT 1 ME . D» CONM OAUX ) 


TC( J)=0. 

M ^R(J)*0. 
IF(DONE) GC 


(DONE) GO TO 998 
INSECTS. F ALSE, 

KeK+ 1 
S(K J=S(J) 

S( J)cO. 

TC( J)=0. 

40 P(J)=0, 

or TO 10 

' ****THIS BPANfF IS USED FOR THE NON-OVErlAP 
COFFlTIpN OF SIGNAL BURSTS**** 

25 IFdHEORYjcnL 1 DE AL ( KH AX* I PO I NT . I T AU, I ©A K P# NT OT* MT OT/ I R 
*DEAD) 

Do 60 Jel.KMAX 

IF C PLOTS) CALL F I OT ( M* J* KmaX. OT 1 ME . 0# CON* I OAUX ) 

Tr ( J)=0. 

IMPDNE) GO TO 996 
InT FRsKMAX 

10 m IMEsINTER+FTIME 
Gp TO 999 

996 T OT IME=IDTV *DT 
PERMf, sNTOT/TOIImE 
tvRl TL ( lOUT* INp) 

iFCPNt) hRlTinrUT* 100)NTOT#MTOT 

lOC FORMATt//” THE NT0T= **120*/* THE MTOT* **I20) 

IF (PNE )S10P 
HR1TL(I0UT*1?0) 

WRITF (I DOT* 110)(MKXfI)*I*I = l*IP) 

110 F DRMAT(5(?X/ J7/**I=*I2)) 

120 F0RHAT(//* ACCUMULATOR SUM FOR I=1*IP IN N( K ) N( K- I ) */ ) 

STOP 

End 

Sub ROUT I KF SIGS( AMP.KMAX/IO, IR1,1R2*D*C0N. AME#F0*KAvE.C0NST>CONS) 
OJMFNSION HC50) 

Logical rave, ci f st* zeros 
Common si looo), tcuooo)»x(47 )* jv 

Common /ado/maxc.gi.gz 

lF( ( .NOT ,2tR0S).AND,WAVE)S( l) = Eyp(-2,*TS**2/D**2)*( 1 . + AmE* 
C0S<6,263165*F t *TS) )*AMP + SC 1 ) 



?fi96* 
2900* 

291 0*t 
2920* 
2930* 
29i»0* 
29f)0* 
2960 * 
2970* 
29B0*C 
2990*C 
3000* 
3010*- 
3020*C 
3030* 
3040* 
3050*C 

lo70*C 
3080*C 
3090* 
3100* 
3110* 
3120* 
3130* 
3140* 
314b* 
3150* 
3160* 
3170* 
3190*t 
3200*C 

3210* 
321b* 
3220* 
3230* 
3240*0 
3250* 
3260* 
3270* 
3280* 
3290* 
3300* 
3310* 
3320* 
3330* 

mu 

3360*C 
3370*0 
3380*0 
3390*0 
3400*C 
3410* 
3420* 
3430* 
3440* 
3450* 
3460* 

IS58: k 
M88: 

3510* 
3520* 
3530* 
3540* 
3550* 
3560* 
3570* 


|F( CONST ) SI 1 )=cr K'S + St 1 ) 

- IMPULSE RESP0^SE OF PMT FROM -TO TO +10 
Do )0 lsi*lVAL 

J=I-I0-1 

H( 1 ) = (1-ABS(F LOATC 1P*DT)*GE 

AsPAND( Jh2# JVj.02) 

On 20 K= 2 fKMAX 

Ts(K-l )*0T + TS 

- IF THIS ISNT A F U-LLRC ZEROS C DM0 1 T I ON ) # AND WE WANT THE 
-classical, WA yl F DRM THEN CALCULATE LAMBDA, 

IF f ( ,NOT,/ERnS), AMO. WAVE) S ( K ) = E XP ( -2 # *T **2/D**2 ) * ( I , + A ME* 
0nSf6,263ie5*Fl *T))*A mP+S(K) 

- IF THERL 15 A CONSTANT BACKGROUND THEN ADD IT INTO LAMBDA 
Ir (CONST )S(K)»fONS + SCK) 

Ss 2=(DT*(SCK)+5( K-1 ) )/2.0)+5Sl 

- CALCULATE THE INTEgF^AL 
SH = (SS1 + SSP/2.C 

- ADD UP THt PHlTONS USING ThE CLASSICAL WAVEFORM 

- THIS IS the RAFDOM PHOTON COUNT (TC) 
i lEt^tGTiSHlGr TO 20 

TC(K-1 UtC(K-1) + 1 
ZsZ + RANd( 1' JV^(>2 ) 

Go TO 1 
^0 SSlcSS2^ 

Do 50 IsWKMAX 

If ( 1R1 ,EQ»4 )G(' TO 5l 

MAXsTC(l) 

TCR( I >*0.0 

IFCMAX,E0,0)G0 TO 50 
Do 60 J=1.F'AX 

60 TCK(1>=RAND(IR1/JV,G1)+TCR(1) 

1 VAL=I0*2+1 

-^CALCULATE THE SIMULATED SIGNAL BURST 
Do 30 L= 1 #KMAX 
OO 40 l=i;iVAL 

Jal-lo-i 

lF( ( L-J) ,LT , 1 , OF , (L-J) , GT. KmAX>gO TO 40 
R(L)=H( I )*TCR(L -J)+R(L) 

40 CONTINUE 
30 CONTINUE 

Return 

ENO 


IS UStD 

- DEPENDING ON THE FIRST INPUT 

Calculate ■ 


FOR 


TO calculate 
'^FXPONEN?! AL 


... AN . ... 

calculate A KAYLIEgH 
calculate a GAUSIAn 
give a constant=i,o 


„R.v, 

R. V. 


I* 1 
= 2 , 

= 3 , 

J=2051*J " ^ 

J=MPD(J, 4194304) 

XsJ/ 4194304,0 
Go TO (1^2^2»4)M 
4 RANDbI.O 
Return 

RAYsRAND 

J=2051*J 

J = MnO( 4194304 ) 

X=J/4194304,0 

RANp=G*(KAY*C0S(6.283l653*X) )/S0RT(2/3« 14159265)+1 ,0 

K t T I .» R N 

end 


A random variable 

R.V. 


14159265) 
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3580* Subroutine^ ideai (inuh,i point# iT«u»iQ,iP»NTOT^MK,iR#DtAD) 

3590*c- THIS SUBROijTINF SIMULATES AN IptAL PROCESSOR 

3600* LOGICAL ZEROS# U SECT*DEDtIM,ONE#DONE 

3610* Common /SINE/MHVCPOT^ONE^DONE^IOTV/IOTTOT^hNTOT , 

3620* Common /idea/tcf(1o6o)#zeros#n(?ooO)#itcr(io)#insect 
3625* Common BULKIOOC )#TcU000)#BUL2(AB) 

3628* Common /adu/maxc #gi#G2 

3630* integer dead 

36A0* NUM*INUM 

3650* In • *0 , 

3660*c-'0N the zeros cONDl T I 0N( SpAcE BETWEEN THE BURSTS) CONFUTE 
3670*C- MTOT with VALUIS OF ZERO 
3680* IF fZEROS)GO TO 9 

TO t . 

If ( INSECT )num=mm+kmAx 
D o 2 J=1#NUH 

JlsNUM-J+1 

2 TC( Jl+IR)sTC( J1 ) 

Do 3 Jsl#IR 

3 TC( J)sITCp( J) 

1F( INSECT )NUM = MiM-KMAX 
1 NUM=NUM+lR 
I=MOD(NUM# ITAU) 

jfyww 1F(I*E6 jU)G0 to 4 
3800*C- IF OIVISIUN RhSULTS IN A REMAINDER (I) THEN SOmE 
3810*C- OF THE RANDOM FHOToN COUNT ARRAY £TC ) HAS TO^ 

3620*C- be SAVED AND TpEN uSEQ ON THE NEX* cALL TO IDEAL 

3830* DO 5 J=1#I 

3840* 5 ITCR( J)sTC(NLP-l+j) 

3850* 4 IR=I 

3660* NkMAXsNUM/ITAU 
3870* lF(NKMAX*EO.O)Ff TURN 
IVAL = NKMAX+IP0UT-1 
Dn 30 KslPyiNT# IVAL 
DO 35 I = 1#ITALI 
IF(DEDTIm)GO to 33 

39?0*C- IF A DEAD TIME CONDITION EXISTS THEN SKIP THE SUM 
3930*C- OF N(k) FOR THIS DT 


.91 
’ 0 ( 
3710* 
3720* 
3730* 
3740* 
3750* 
3760* 
3770* 
3780* 
3790* 


3680* 
3690* 
3900* 
3910* 
19 
19 


3940* lPHOTSeTC(l+(K-IPOlNT)*lTAU) 

3950* IF( TPHOTS.GT.MAXOIpHOTSsMAxC 

3960*C- MAXC PHOTONS IS THE LIMIT OF THE COUNTER IN THE 
3970*C- IDEAL PROCESSIF 

3980* , IF(TPH0TS.GE,1, AND. dead, GT,0)0EDTIM«. TRUE, 

3990*1- IF A DEAD TIM^ FROM SATURATION IS t 6 BE f VALUATED# AND The 
4000*C- threshold of 1 PHOTON IS CrOSSFD THEN SET THE LOGICAL 

4oio*c- dedtime condition and Skip dead numbef< of dt s 

J8?8; 


8818 : 

4060* 

4070* 

4080* 


TO 


35 


DEDTIMs, FALSE. 

35 CONT!nUE 
4090* 30 continue 
4100* Go TO 20 

9 IF(IR,EQ.0)G0 to 10 

num=num+ir 
D o 12 J=1#IR 

12 N(IPOINT)=N(IPOInT)+ItCR(J) 

10 IF(NUM.GE,1P)G0 to 14 
Dn 11 J=IR+1#NLF 

11 TTCR(J)=0 
IR=NUM 
N( IPOINT)sO 

return 

14 I-M00(NUM# IT AU) 

IF(I.EO,0)6o TCi 16 
DO 13 Jsl#l 
j3 ITCR(J)*0 


4110* 

4120* 

4130* 

4140* 

4150* 

4160* 

4170* 

4180* 

4190* 

4200* 

4210* 

4220* 

4230* 

4240* 

A250* 

4260* 
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NkmIxsJum/itau 


!vAL*l4K!^AX+iiJQlNT-l 

i? 

20 IP01NT»IPQINT+NKMAX 
OU <tO L»l^lPOINr 

IF(ONE)GO ro 43 
00 42 J»IpIF 

42 MKX( J)«MKX( J)+N(L)*N(L-J) 

43 lr((L-IO),LT,l)<jO TO 44 
f<KsNa,)*CN(L-lP)-NCL-»ia) )'fMK 

44 IDTV1»L*ITAU 

iF(ioTvH'ioi\/.GE.ioTronGO ro v? 

NTOT»MTOT+MCL) 
iFCNTOT.GE.MNTQT )G0 IvJ V9 
40 CONTIi^UE 
00 4& K»1»IP 
KIslPQI(4T-iP + K 
4b i4(<)sN(Kl) 

IPOlNTsIP+l 
00 30 K*IP01NT^KI 

toTVafJfv+lOTVl 

5y^Oo3E*.TRUE. 
loTVslorV+IOTVl 
RETURN 

: {3^ IlgSit M 

• AND THE random PHUTON COUNT I*) 

883383 'simmmu 0)»X(4T),JV^PHI,PL0 
MsM + 1 
SVALsS( J) 

RVAL*R( J) 

if(Sval.lt.plo)sval=plo 
lFCRVAL.aT.PHl )R vAL*PHI 

iffMlsJaf ?4?)5sSs4/ 

JMAXsHSS 

lF(M^jGT,MSS)JMAx*MS 


(TV+IOTVI 
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‘»6ao* 

j6;>o* 

4660* 

46/0* 

4660* 

4gy0* 

4/00* 

4/10* 

4720* 

4760* 

4740* 

4760* 

4760* 

4770* 

4760* 

47V0* 

4300* 

4310* 

4820* 

4860* 

4840* 

4860* 

4860*C 

4870* 

4860* 

48V0* 

4900* 

4910* 

4920* 

4960* 

4940* 

4960* 

4960* 

4970* 


LUrtl*MC*l 
KHI CiHsJMAX-1 

it 

Ju'yS'SliJiJg'’ ” 

2S XEA1)=1H» 



■JHlls“sU8RJUTINFpLUrs''s?ACE'(AfLL;S^ ^ 

3SM + 1 bKflOt CMLLEKj IM 3ETwEEN 3URS>IS 

If (M,£ o, 5 jau TO 23 


4Ki rE(IA0X»65) 
NR,HAT(1X»« 
au ro 70 


lio ro 70 w*ooo 

3sO 0.000»»,3XxL 10, Ax »* + #«) 


‘♦3 

MsO 

7o^ketok^ 
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APPENDIX C 


VARIANCE OF THE DUAL CORRELATE 
AND SUBTRACT ESTIMATOR 

In this appendix we wish to determine a formula for the 

mean-square deviation of the estimate ill defined in the text 

equation (56) at value of delay (near 11^/^ where the expected 

value of M is zero. It was our hope to do this for the 

pq 

general signal model which we have presented. To complete such 
a task requires the use of fourth-order moments of the classi- 
cal signal process X(t), and it was for this reason that the 
derivation in the last part of Appendix A was undertaken. We 
did not have time during the contract period to evaluate and 
use the fourth order moments of the general signal. For this 
reason we have restricted this analysis to the case where 
steady background light is the predominant source of variability 
error. (Steady light adds variability error even though it 
cancels in the mean. ) This simplifies the problem because of 
the simplicity of the fourth order moment equations for steady 
light (homogeneous Poisson counting process). 

We define here as* 
pq 


V = ^ ”k (Cl) 

"k ‘ “k“k+p - “k“k+q 

where the summations will all be from 1 to N unless otherwise 
noted. We have 


<m^> = <nj.xn^^p> - <ni^><i>t+q> 


( 02 ) 


‘^This derivation was performed with plus signs in the 
delay subscripts instead of minus signs. There is no differ- 
ence in the results. 
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where n = XAt = nArP/hv. From (C2) we have 


<M„^> = 0 

pq 


(C3) 


we wish to evaluate var -■ 

pq pq pq pq 


<m: 


pq 1 j ^ X X j 


(C4) 


In (C4) the product of sums was expanded and the order of 
expectation and summation interchanged. The i = j terms are 
separated because they behave differently. 

The terms in the first summation give 

<ml > =. (05) 

= <"?>!<n|+p> + <nf^.q> - 2<n^^pxn.^^>] 

= (n + n^)(2(n + n^) - 2n^) = 2n^ +-2n^ 


where the theory of homogeneous Poisson processes has allowed 
the factorization due to independent" of nonoverlapping count 
intervals* and from which we know that 

<n^> = n + n^ (C6) 

Next we must evaluate the i ^ j terms in (C4). There 

2 

are N - N such terms but many of them are zero. We have 


* 

2 At this point for the general nonsteady signal evaluate 
<m^ > with X(t) conditionally given and then evaluate fourth 
oraer moments of the process X(t) which is also Poisson. 
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<m.m .> 
1 3 


(C7) 


= <n.n.^pn.n.^p> - <n .n.^^n 


-<n.n.. n.n.. > + <n.n.. n.n.,„> 

1 i+p 3 j+q 1 i+q 3 j+q 


When all the subscripts are unequal we obtain 


— 4 — 4 — 4 — 4 « 

<m.m.> = n- n - n +n = 0 
1 3 


We are restricted to i j and p ^ q 7 ^ 0. If we examine the 
matrix of products m.m., the only allowable products for which 

J- 

some of the subscripts of n in (C7) are equal, excluding the 
i=j case, are found on diagonals parallel to the i=J diagonal. 
These diagonals are i=j ±q, i=j ±P, and i = 3 ± (p-q). 

We examine these diagonals separately: for i = 3 + q 


<m . m . > 
1 3 


<n .. n. . n.n.. > — <n., n..^ n.n., > 

j+q j+q+p 3 j+P j+q o+2q 3 3+p 

2 2 
- <n.. n.. . n.> + <n., n.,o a.> 

3 +q 3+q+P 3 J+q 0+2q J 

— 4 — 4 - 2.-2 >— - 2.-2 

n - n “ (n + n )n + (n + n )n =0 


(C 8 ) 


Similarly for i=j -q, i=3 ±P> the result is zero by sub- 
traction. For the diagonals i = 3 * ± (p-q), however, we obtain 
after substitution and evaluation 


<m.m.> = -n^ (C9) 

X J 

For large N we may neglect end effects and observe that there 
are approximately 2N terms which result from these two diagonals. 
We may now evaluate (C4) using (C5) and (C9) as 

<M^ > = N(2)(n^ + n^) - 2Nn^ (CIO) 

= 2Nn^ = 2N(XAx)^ * Var(M ) 

ir 

This is the desired result . 
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APPENDIX D 


SYSTEM DESIGN OF AN ADVANCED PHOTON COUNTING 
PROCESSOR FOR LOW-LEVEL LDV SIGNALS 

Introduction 

A problem confronting us in the present design was one of 
complexity and cost. A two channel duplication of the AEDC 
processor [1] with the extra requirements of dual-channel 
advanced-concept operation would have been prohibitively 
expensive. To solve this problem we have made several sig- 
nificant changes, some of which utilize the power of a high- 
speed computer controller. 

Design Approach for the Advanced Processor 

Major cost savings are associated with exact duplication 
of circuit layout. For this reason, the dual channel system 
is designed so that it may be operated either as two com- 
pletely separate identical units in one rack; or as two inde- 
pendent units with synchronized data (for cross correlations) 
or as a one channel system with a synchronous normalizing chan- 
nel (same system clock and n^^ sequence but separate delays). 

The counter/timer (C/T) functions needed for system control 
were obtained in the AEDC system [1] from a $1500 laboratory 
counter, selected because of the availability of options ($285) 
for computer control and read. It is considerably less expen- 
sive to include the C/T functions in the special-purpose hard- 
ware to avoid both the cost of two units ($3570) and the asso- 
ciated computer cables, I/O cards, etc. This puts an even more 
stringent requirement on reducing complexity of other circuits 
to give more room on the wire wrap panel (162-180 IC sockets 
unless multiple panels are used). The following system concepts 
have been incorporated. 
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1. Eliminate all front panel controls and displays. 

The system thus utilizes the display and command 
of the controlling mini-computer and cannot func- 
tion separately. This saves LED driver integrated 
circuits as well as the front panel itself. 

2. Use a newer, lower cost voltage variable oscillator 
for the system clock. 

3. Replace the laboratory counter with ECL counter/ 
timer (C/T) circuits. 

4. Use an external, separate heavy duty ECL power supply 
(=: 16 Amps) to save space in the system enclosure and 
avoid a large heat source in the -enclosure. 

5. Reduce the maximum count to 3 bits instead of 4 
(see following justification). This simplifies 
all of the circuits and reduces package count. 

6. Reduce the number of control and read data circuits. 
(See following discussion. ) 

7. Do not require the 3 bit counter to be selectively 
saturable. Let it saturate at 7 (111). 

8. Remove the single-clipper circuit and the associated 
multiplexer. 

9. Limit the accumulator to 15 bits + sign to be com- 
patible with the 16 bit word of the computer. Read 
the accumulator often enough to avoid overflow. 

10. Replace the 4 bit subtracter with a 3 bit subtracter 
adder. The add function could be used in one of the 
normalization schemes. 

11. Include an experimental analog servo loop for zero 
adjusting the system clock to the proper multiple 
of the mean signal frequency. 


Reduction of Counter Bit Number 

The maximum periodic rate of discriminator output pulses 
is 120 MHz (dead time > 8nsec). The largest random rate is 
usually less, say 70 MHz. In order for the count of 7 to be 
executed therefore, we assume the count interval to be larger 
than 0.1 ysec (1/fQ > 0*4 or 0.8 ysec depending on the delay 
choices). Thus we may have the count exceeded for signal fre- 
quency f^ <2.5 MHz or 1.25 MHz. We observe, however, that a 
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count of 7 in a quarter of a cycle is a peak rate of 14 elec- 
trons/cycle. This is where the burst counter can become 
useful . Thus a 3 bit counter is adequate because when larger 
signals are present a burst counter processor should be used. 

If the larger signals are unwanted because they occur from 
larger particles, then saturating the count at 7 will reduce 
their effect in favor of smaller signals. 

Read logic for photon statistics .- The 16 bit word length 

and the high speed of the H.P. 2100 computer will allow straight 

2 

forward READ circuitry. Each channel will have a 16 bit T L 
latch connected to the 15 bit plus sign acctimulator . In order 
that the counter/timer and the multiple interval statistics 
may be read, selectable steering gates will be used so that the 
first 6 bits can be connected to the output of the multiplier. 
Sequential numbers may then be read as follows: 

1. Set maximum delay in A+B and set one side of multi- 
plier to 1. 

2. Run data and clock. 

3. Stop clock. 

4. Read multiplier out with commanded single clock 
advance. 

This will produce a string of more than 20 sequential numbers. 

In order to read the C/T total (8 digits) 32 bits is required 
(BCD code). This will be accomplished by sequentially reading 
two 16 bit words time-multiplexed to the TTL output port. 

Description of Circuits 

This section describes the subsystem circuits and compo- 
nents which comprise the dual channel photon processing system 
one channel of which is shown in Figure D1 . 

External subsystems .- The external subsystems are the 
following: 

1. Two photomultiplier tubes and associated housings 
and power supplies. 
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2. Six precision 50 connecting cables. 

s 

3. Two preamps, 

4. Two amplifier/discriminators. 

5. One NIM bin (rack mounted). 

6. Two voltage controllable oscillators (1-200 MHz) 

7. Two rack-mounted integrated circuit power supplies 
and connecting cables. 

8. A computer with two microcircuit I/O cards (16 bits 
in and 16 bits out for each with device command line 
and device-ready flag line), and two 36 pair twisted 
lead cables with connectors. 

Three bit counter . - This is a dual section counter which 
stops and holds the count of seven instead of turning over to 
zero and continuing. Two counter sections alternate so that one 
may have data transferred and be reset while the other one is 
counting. In order to avoid the possibility of counting a 
border line event twice, a dead time between count intervals 
approximately 2 nsec will be incorporated. The dead time 
between input pulses (from the discriminator) will be ^ 10 nsec 
with a design objective of 8 nsec. The alternate count inter- 
vals will be equal (design objective). A control bit allows 
the counter output to the delays to be set to zero (data gate). 

Delay sect ions . - The two delay sections A and B are speci- 
fied as A = {0,1,2,3,4,5,6,8,9,10,11,12} and B = {0,2, 4, 8}. 

This set allows sequential autocorrelations with delays up to 
20 in addition to the dual correlate and subtract and the 
normalizing modes. 

For dual channel operation with one channel normalizing, 
a multiplexer is provided which allows the data from the other 
channel counter to be selected as delay input. Since the sep- 
arate channels are identical physically, each has an input and 
an output to the other channel. 
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Four control bits select the A delay. Two control bits 
select the B delay. Ah additional bit selects the output of 
the other data channel instead of the counter output of this 
channel . 

Adder/subtractor & mode selector . - In the AEDC unit an 
option for subtracting a constant from each nj^ was included. 

The purpose was to avoid overflowing the accxamulator in auto- 
correlate mode when a signal with large mean is encountered. 

The ability of the computer for high-speed read to DMA 
avoids the necessity of subtracting a constant, since the accu- 
mulator may be read often to avoid overflow. The two control 
bits allow selection of add, subtract, or add zero. The adder 
subtracter output is four magnitude bits plus sign. 

Multiplier .- The multiplier has seven bits plus sign out. 

A control bit allows the undelayed path input to be set to 001 
instead of nj^. 

Accumulator . - The accumulator has 15 bits plus a sign 
bit. It can be reset by a single pulse. It will be implemented 
in 2's complement; the computer software will convert to sign 
and magnitude. The reset pulse must also transfer the accumu- 
lator values to latches in the output port to the computer. 

Counter/ timer . - This section replaces the external counter/ 
timer used in the AEDC System. It consists of two 8-decade BCD 
counters and associated input selectors and controls. Each 
counter may select as input either the precision 1 MHz oscil- 
lator, the system clock, or the second discriminator output 
(external input). A fast prescaler selects divide by 1, 2, or 
5 for counter #1. Counter #2 is multiplexed- to the 16 bit out- 
put port (1st 4 significant decades and 2nd 4 significant 
decades separately selectable. ) The first counter produces 
output pulses at decade counter #1 selector. The package count 
does not include multiplexing and read outputs for counter # 2 . 
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A gate generator block does the following: It produces 

the pulse which is applied to the accumulator and the computer 
flag; it programably reads and resets the C/T; and it produces 
a gate pulse between two successive counter #1 output pulses 
for application to the counter #2 control gates (counter #2 
stops and holds when gate goes back down). The programable 
states include: a) read and reset counter #2 at the entry of 

the control word (CW) and at each output pulse from counter #1 
(PFC); b) read and reset counter #2 CW but not PFC; c) read 
and reset counter #2 at PFC but not at CW; d) do not read and 
reset counter #2 at either PFC or CW. (This allows, by software, 
for counting the photon rate over the duration of an autocorre- 
lation sequence and then stopping at the end of a number of 
smaller intervals. ) 

The control bits needed for the C/T are as follows: 

Three bits are used to select which of 3 inputs goes to each 
of the two decade counters. (Assumes both will not have same 
input.) Three bits will select the decade output of the decade 
scaler (counter #1). Two more bits select the prescale divisor 
(1, 2, 5). Three bits have been allowed for the four read/ 
reset states (2 bits) and an extra control bit (spare). Finally, 
there is one bit which allows application of the control gates 
from the other C/T to the accumulator read/reset line instead 
of this channel's C/T. (This is for cross correlation with 
synchronous accumulator read/reset.) 

System clock .- The system clock accepts a periodic wave 
form with 1 positive-going transition per cycle) and shapes this 
into a periodic pulse train at the ECL voltage levels. It 
includes buffer gates for proper fan out (5 packages). It 
includes controls which allow the clock to be stopped cleanly 
and a one-shot clock pulse generator which can be activated 
by computer instruction. This feature allows sequential n^^ 
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values to be read out of the system for multiple interval 
statistics. The control bits are one for on/off and one for 
single pulse. There is also an ungated clock line which is 
used as an input to the C/T. There is a control bit which 
allows the other system clock to be selected (when the other 
channel data is used, for example). 

Output port . - Each channel has one 16 bit TTL compatible 
output port and a 1 bit flag pulse line. The output port 
includes a multiplexer (two 16 bit sections of the C/T, the 16 
bit accumulator output and a 7 bit output from the multiplier), 
a 16 bit latch, and ECL to TTL voltage level translators. 

The output multiplexer requires 2 bits of control to select 
one of 4 outputs. 

At the present time it is not clear which of the following 
approaches could be utilized: a) all ECL construction with 

only 17 bits of ECL/TTL translator as the output; b) a separate 
panel section of TTL circuitry which includes most of the C/T 
and the output latches for the output port. The b) approach 
would require less expensive IC's and less power for part of 
the system. However, the cost of a separate panel and the 
panel interconnections may make a) preferable. 

Oscillator control .- This subsystem provides automatic 
fine tuning of the variable system clock for the determination 
of the mean velocity. The mean value of the multiplier output 
is negative while the clock frequency is too high. Thus the 
intent here is to use fast digital/analog conversion and 
analog integration with controllable reset, integration, and 
hold states to provide a control voltage to the external volt- 
age controllable oscillator. Two control bits are required. 

This feature is a research item. It may later prove more 
advisable to use the D/A converter with a portion of the 
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computer control word bits so that the mean velocity selection 
becomes a software task. 

Input port .- All machine control is accomplished by 30 
bits of information and the time at which certain of the bits 
are changed. These bits are latched into the processor input 
port 15 at a time by a single line command pulse from the com- 
puter. All time-critical bits of control are included in con- 
trol word 0. (One of the sixteen bits from the computer is 
the address of the control word 0 or 1.) Table D1 provides a 
tentative assignment of control bits. The bit number refers 
to the power of 2 in standard binary format. 

Package count . - An integrated circuit package count esti- 
mate of 147 IC's made for one of the two identical channels. 

The estimate assumes the use of a panel with 162 sockets and 
design for each subcircuit with the same approach previously 
developed for the AEDC unit . This leaves 15 spare sockets for 
flexibility in design and/or additions. A slightly improved 
approach has also been identified which uses 2 bit arithmetic 
logic units (ALU) instead of 4 bit ALU's to reduce circuit 
speed limitations. This approach would actually produce more 
useable sockets because the panel would have 180 standard 
sockets instead of 150 plus 12 4 bit ALU sockets. 
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Address Zero 


Bit Function 


Bit , Function 


0 Clock Gate 

1 Clock One Shot 

2 3-Bit Counter Gate 

I! 

g j A Delay 

6 

g [ B Delay 

9 Multiplier 

10 > 

11 > Output Port Multiplex 
12 ) 

13 [ Accumulator Reset Control 

14 ) 


0 ) 

1 [ C/T Input Selection 

2 ) 

3 ) 

4 [ C/T Decade Select for Counter #1 

5 ) 

® I C/T Prescale Select 

g [ Adder/Subtractor 

10 Accumulator Control from Cross Channel 

11 n from Cross Channel 

12 C^ock from Cross Channel 

13 » 

f Oscillator Control State 


TABLE Dl. CONTROL WORD FORMAT 
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